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I.  INTRODUCTION 

DTIC  QUALITY  INSPECTED  2 

The  performance  of  rotary-wing  aircraft  is  limited  by  a  host  of  aerodynamic  problems. 
Foremost  among  these  problems  are  transonic  flow  (which  is  the  primary  drag  source  on  high¬ 
speed  propellers  and  advancing  helicopter  rotor  blades)  and  stall  (on  retreating  helicopter 
blades);  but  there  are  many  other  flow  problems — especially  on  helicopters,  which  operate 
under  a  wide  variety  of  conditions  and  constraints.  For  example,  a  helicopter  often  operates 
under  conditions  of  flow  separation  and  in  its  own  wake.  Because  of  such  flow  complexity 
the  aerodynamic  design  of  helicopters  is  traditionally  an  empirical  craft  that  often  relies  more 
on  experience  and  test  than  on  detailed  analysis.  This  approach  requires  a  high  degree  of 
design  conservatism.  In  spite  of  a  cautious  design  philosophy,  rotorcraft  development  often 
encounters  unexpected  and  sometimes  dangerous  aerodynamic  problems.  Such  occurrences 
will  persist  as  long  as  we  require  new  designs  but  lack  the  ability  to  predict  and  control  the 
details  of  rotorcraft  flows. 


The  ability  to  design  rotorcraft  with  confidence  requires  a  new  order  of  aerodynamic  pre¬ 
dictive  technology  that  is  both  true  to  the  basic  flow  physics  and  readily  usable  by  industry. 
The  burgeoning  field  of  computational  fluid  dynamics  (CFD)  holds  the  promise  of  provid¬ 
ing  the  necessary  new  tools.  In  the  fixed-wing  community,  CFD  has  indeed  revolutionized 
aerodynamics — even  to  the  extent  of  permitting  the  freezing  of  important  design  features  on 
the  basis  of  computations  alone.  The  rotorcraft  community  is  currently  far  from  this  level 
of  confidence,  because  the  basic  flow  phenomena  are  more  complex,  diverse,  and  interdepen¬ 
dent  (i.e.,  many  phenomena  occur  simultaneously  and  affect  each  other  directly  or  indirectly 
through  their  effects  on  the  flight  state). 

Some  of  the  basic  flow  phenomena  are  illustrated  in  Fig.  1.  The  first  area  of  aerodynamic 
concern  on  a  helicopter  is  the  design  of  the  main  rotor,  to  maximize  hover  and  forward-flight 
performance  and  minimize  vibratory  loads  and  noise.  These  are  conflicting  requirements, 
which  necessitate  difficult  design  compromises.  For  example,  the  choice  of  rotor  tip  speed  is  a 
compromise  between  minimizing  transonic  flow  effects  on  the  advancing  rotor  and  stall  effects 
on  the  retreating  side.  However,  this  compromise  is  strongly  affected  by  the  choice  of  airfoils, 
planform,  and  twist,  and  the  propulsive  requirement.  The  determination  of  these  features,  in 
turn,  depends  on  the  fuselage  drag  and  the  rotor  flow  environment,  which  is  strongly  affected 
by  the  rotor  motion  (flapping  and  elastic  deformation)  and  its  wake  structure.  Therefore, 
a  rotor  analysis  is  a  holistic  process  that  involves  the  blade  aerodynamics,  dynamics,  elastic 
properties,  and  fuselage  effects.  The  analysis  of  various  flow  features  are  central  to  this  process. 
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Rotor  CFD  applications  have  not  yet  reached  the  level  of  sophistication  that  characterizes 
fixed-wing  computations,  which  include  fairly  extensive  Reynolds-averaged,  Navier -Stokes 
solutions.  However,  these  computations,  which  must  in  general  be  steady  because  of  machine 
limitations,  suffice  for  most  fixed-wing  applications.  Comparable  viscous  computations  are 
prohibitively  costly  for  helicopter  rotor  analyses,  because  these  must  usually  be  unsteady. 
Furthermore,  there  exist  more  immediate  (and  non-viscous-dominated)  computational  issues 
that  remain  to  be  solved.  These  include; 

1.  ROTOR  WAKE  PREDICTION.  The  prediction  of  the  rotor  wake  structure  (its 
geometry  and  circulation  distribution)  is  one  of  the  most  important  rotor  flow  issues. 

For  hover,  the  wake  is  the  single  most  important  flow  issue,  because  it  is  the  primary 
determinant  of  induced  power.  However,  the  ability  of  any  code  to  predict  this  wake  with 
sufficient  accuracy  is  a  subject  of  great  debate.  Certainly  there  exists  no  industry  code 
which  can  reliably  predict  the  detailed  effect  of  blade  geometry  (planform,  for  instance)  on 
the  wake.  This  is  because  almost  all  current  analyses  are  boundary-integral  codes  which 
represent  the  rotor  as  a  lifting  line  and  do  not  treat  the  3-D,  non-linear  compressible 
blade  flow. 

For  high-speed  forward  flight,  the  wake  is  somewhat  less  important,  because  it  convects 
more  rapidly  away  from  the  blade.  However,  the  wake  is  still  close  enough  that  a  detailed 
knowledge  of  its  nature  is  required  to  predict  vibratory  loading.  Also,  it  is  generally 
conceded  that  we  do  not  understand  the  structure  of  advancing  rotor  wakes  (because  of 
both  experimental  and  computational  difficulties)  and  we  cannot  reliably  predict  rotor 
vibrations.  For  lower  advance  ratios  and  for  descent  conditions,  the  rotor  and  wake 
closely  or  directly  impinge  on  one  another.  No  general  analysis  of  this  condition,  called 
Blade/Vortex  Interaction  (BVI),  exists. 

2.  COMPRESSIBLE  AERODYNAMICS.  At  high  advance  ratios  the  predominant 
rotor  drag  source  is  due  to  transonic  flow  at  the  tips  of  the  advancing  blade.  On  a 
high-speed  tilt-rotor  this  transonic  flow  could  conceivably  involve  the  entire  rotor.  The 
high  drag  results  from  the  occurrence  of  supersonic  flow  and  shocks.  That  such  flows 
are  intrinsically  unsteady  and  strongly  three-dimensional  demonstrates  a  weakness  of  the 
present  use  of  2-D  measured  airfoil  data  to  predict  rotor  loads  and  performance. 

To  date,  the  most  conspicuous  strides  in  rotor  CFD  have  been  in  the  development 
of  methods  to  solve  the  compressible  flow  problem  using  potential  techniques.  Although 
some  standardized  codes  have  now  been  developed,  numerous  problems  remain,  including 
the  development  of  suitable  methods  to  treat  the  boundary  layer  and  coupling  methods 
to  iiitcrface  these  solutions  to  a  global  rotor  analysis. 

3.  INTERACTION  PROBLEMS.  Rotorcraft  are  unique  in  the  large  number  of  compo¬ 
nent  and  flow  interactions  that  occur.  The  main  rotor  can  interact  closely  with  its  own 
shed  wake  (a  source  of  noise  and  vibration),  whose  location  is  affected  by  the  fuselage 
flow  field.  The  fuselage  upwash  also  introduces  rotor  vibratory  forces.  In  addition,  the 
main  rotor  and  tail  rotor  can  interact  in  liover  to  reduce  the  main  rotor  efficiency. 

Most  of  the  above-mentioned  problems  do  not  involve  strong  viscous  forces  and  can  be 
treated  by  potential-based  methods.  Nevertheless,  unsteadiness  and  structural  and  dynamic 
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coupling  makes  the  helicopter  problem  a  very  large  one.  Moreover,  the  helicopter  problem 
does  involve  very  significant  viscous  complications  which  require  resolution.  These  include 

1.  RETREATING  BLADE  DYNAMIC  STALL.  Blade  flapping  and  cyclic  input  (re¬ 
quired  to  satisfy  propulsive  requirements)  drive  the  retreating  blade  to  very  high  lift  co¬ 
efficients.  Vortex  interactions  on  the  retreating  side  also  exacerbate  the  problem.  When 
stall  does  occur,  it  happens  after  a  dynamic  delay  and  is  far  more  violent  (inducing  very 
high  drags,  pitching  moments,  and  flutter)  than  steady  stall  is.  This  defines  a  hard  perfor¬ 
mance  limit.  The  process  is  also  strongly  3-D,  but  no  method  has  yet  been  demonstrated 
which  consistently  predicts  even  2-D  unsteady  stall  behavior.  Development  of  the  nec¬ 
essary  viscous  flow  methods  is  an  important  long-term  goal.  In  the  interim,  however, 
empirical  engineering  stall  methods  will  be  required  and  it  will  be  necessary  to  find  ways 
to  include  these  in  our  simpler  (potential-based)  CFD  flow  analyses. 

2.  FUSELAGE  FLOWS.  At  high  speeds,  the  largest  helicopter  power  sink  is  the  parasite 
drag  of  the  fuselage.  Indeed,  alleviating  this  drag  (and  hence  the  propulsive  requirement) 
is  a  most  potent  way  of  minimizing  the  retreating  blade  stall  on  the  rotor.  The  devel¬ 
opment  of  methods  to  predict  the  details  of  the  fuselage  flow  therefore  assumes  great 
importance.  Helicopter  fuselages  are  often  unavoidably  bluff  and  this  causes  separated 
flows.  To  date,  our  ability  to  predict  the  drag  of  such  bodies  is  limited  and  requires 
empirical  input.  This  ability  is  another  long-term  flow-prediction  goal  and  one  that  can 
benefit  from  fixed-wing  aircraft  experience. 

Thus  a  wide  variety  of  near-  and  long-term  rotorcraft  flow  problems  await  resolution  by 
CFD  methods.  There  are  also  many  CFD  rotor  applications  in  current  use.  This  paper  will 
introduce  prospective  users  or  developers  to  some  of  these  methods  and  to  the  issues  involved 
in  their  use. 

This  exposition  will  begin  with  a  description  of  the  various  flow  equations  and  the  physical 
approximations  they  embody.  The  basic  methods  of  solving  these  equations  will  be  discussed. 
This  discussion  will  be  aimed  at  methods  that  can  enter  common  engineering  use,  now  or  in 
the  near  future,  and  therefore,  will  be  mainly  confined  to  the  various  potential  approximations. 
Many  methods  have  already  been  developed  to  treat  various  aspects  of  rotor  flows,  and  some 
of  these  are  in  wide  use.  The  discussion  of  various  potential  methods  will  include  the  means 
to  combine  these  different  analyses  to  obtain  appropriate  hybrid  codes.  A  novel  approach 
to  obtaining  a  unified  solution  (rather  than  a  hybrid)  will  be  discussed  in  the  context  of 
performing  hover  computations.  Methods  will  also  be  discussed  for  coupling  CFD  analyses 
to  comprehensive  codes  to  obtain  trim  solutions.  Comparisons  will  be  made  with  data  from 
various  flight  and  special-purpose  experiments  to  demonstrate  the  efficacy  of  the  varying  flow 
treatments.  Finally,  some  advanced  flow  topics  (including  Navier-Stokes  solutions)  will  be 
outlined,  in  order  to  demonstrate  the  future  possibilities  as  newer  CFD  methods  mature  and 
become  practical. 


3 


II.  FLOW  EQUATIONS 


The  complexity  of  rotor  behavior  has  required  the  use  of  a  hierarchy  of  flow  models  that 
are  fast  and  pliable,  and  enable  a  comprehension  of  the  entire  system.  The  many  limitations 
of  these  approximate  models  have  driven  the  development  of  CFD  for  rotor  applications. 
However,  even  the  most  elaborate  and  expensive  of  the  CFD  models  involve  approximations, 
which  break  down  under  limiting,  but  not  unusual,  conditions.  An  awareness  of  the  limits  and 
costs  of  the  flow  models  is  therefore  useful.  The  following  is  accordingly  concerned  with  those 
flow  models  (from  most  comprehensive  to  the  simplest)  on  which  the  new  rotor  flow  methods 
are  based. 


NAVIER-STOKES  EQUATIONS 

Many  of  the  critical  phenomena  of  concern  involve  dissipative  processes  (mainly  viscous 
forces).  These  processes  include  boundary  layers,  separation,  and  vortex  formation.  There  are 
other  flow  processes  that  require  viscous  dissipation  but  are  effectively  described  by  inviscid 
models.  These  processes  include  shocks  and  the  airfoil  Kutta  condition.  Viscous  phenomena 
are  so  ubiquitous,  however,  that  the  need  for  a  general  viscous  flow  model  will  always  be  with 
us. 

The  Navier-Stokes  equations  are  considered  to  be  valid  for  modeling  all  continuum  flow 
processes.  For  a  perfect  gas,  these  are 


mass  conservation: 


momentum  conservation: 


^  d{puj) 
dt  dx 
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djpuj)  djpUjUj)  ^  dTjj 

dt  dxj  dxi  dxj 

energy  conservation: 

d{ph)  d{phuj)  _  dp  dp  dui 

dt  dxj  dt  ^  dxj  ^  dxj 


dqj 

dxj 


(1) 

(2) 

(3) 


where  p  is  the  density,  u  the  velocity  vector,  p  the  pressure,  h  the  specific  enthalpy,  the 
stress  tensor,  and  g  the  heat  flux  vector.  The  thermodynamic  quantities  are  related  by  the 
state  equations 


p  =  pRr  (4) 

and 

h  =  CpT  (5) 

where  R  is  the  gas  constant  and  Cp  is  the  specific  heat  at  constant  pressure.  The  only 
assumptions  to  this  point  are  that  the  fluid  is  a  continuum  and  a  perfect  gas.  (The  latter 
assumption  would  not  hold  under  condensation  or  icing  conditions.)  The  heat  flux  is  given  by 

fVT* 
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where  k  is  the  thermal  conductivity.  The  relation  between  stress  and  rate  of  strain  is  given  by 


where  6ij  is  the  Kronecker  delta,  A  =  — 2/3/Lt  the  bulk  viscosity  (Stokes  hypothesis),  and  /x  the 
viscosity.  This  expression  is  based  on  the  assumption  that  the  fluid  is  isotropic  (stress/rate- 
of-strain  relation  is  not  direction  dependent)  and  Newtonian  (stress/rate-of  strain  relation  is 
linear).  Although  Eqs.  (1-7)  are  believed  to  be  correct,  we  cannot  solve  them  without  making 
major  simplifications. 

Many  numerical  solutions  of  the  actual  time-dependent  Navier-Stokes  equations  have 
been  obtained.  However,  these  are  usually  only  valid  for  laminar  flows.  For  Reynolds  numbers 
greater  than  about  1000,  flow  becomes  turbulent  and  dissipation  becomes  dominated  by  mixing 
(rather  than  molecular)  processes.  The  numerical  resolution  of  the  various  characteristic 
mixing  scale  lengths  requires  an  extremely  dense  grid  if  the  above  equations  are  to  be  used.  At 
the  present  time  such  computations  require  tens  to  hundreds  of  hours  on  the  fastest  machines 
(e.g.,  the  CRAY-2),  for  flows  in  regions  which  are  too  small  to  be  of  use  to  the  aircraft 
engineer.  Clearly,  the  importance  of  the  Navier-Stokes  equations  is  that  they  are  required  for 
the  eventual  understanding  of  turbulent  processes,  and  they  will  motivate  the  development 
of  supercomputers  for  many  years  to  come.  However,  their  engineering  importance  lies  far  in 
the  future. 

In  order  to  derive  a  unified  set  of  viscous  equations  with  application  potential,  resort  is 
made  to  Reynolds  averaging.  This  has  the  effect  of  multiplying  the  number  of  dissipation 
terms  and  introducing  new  effective  viscous  and  conductive  terms  that  require  an  assumed 
mathematical  model  and  defining  relations.  Most  such  mathematical  models  consist  of  a 
simple  algebraic  or  mixing-length,  effective-viscosity  expression.  Such  models  are  probably  of 
qualitative  interest  only  for  the  modeling  of  separation  effects.  However,  these  are  currently  the 
most-used  viscous  models  because  they  are  the  simplest.  For  modeling  extensive  separation, 
more  complex  models  are  required — k  -  e  models,  for  instance.  However,  such  viscous  models 
require  the  addition  of  two  or  more  convection  equations  to  the  total  system  and  reduce 
the  size  of  the  problem  that  can  be  modeled  with  present  resources.  One  frequently  used 
simplification  involves  the  elimination  of  the  streamwise  viscous  stress  terms.  The  resulting 
equations  are  referred  to  as  the  thin-layer  Navier-Stokes  (TNS)  equations.  The  rationale  for 
this  simplification  is  that  streamwise  viscous  stress  terms  are  usually  much  smaller  than  the 
transverse  stress  terms.  Furthermore,  memory  limitations  usually  do  not  permit  sufficient  grid 
points  to  resolve  the  streamwise  viscous  stress  (for  a  3-D  problem).  Therefore,  for  problems 
of  real  engineering  interest,  most  unified  viscous  computations  use  the  TNS  equations  with  an 
algebraic  viscous  model.  Furthermore,  most  of  these  computations  are  performed  for  steady 
flows. 

EULER  EQUATIONS 

The  above  discussion  should  give  some  idea  of  the  magnitude  of  the  simplification  which 
results  when  an  inviscid,  nonconducting  flow  model  can  be  used.  When  a  grid  can  be  sized 
to  resolve  body  geometry  rather  than  viscous  scales,  the  grid  reduction  can  easily  exceed  a 
factor  of  ten.  The  resulting  equations  are 
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mass  conservation: 

dp  dipuj)  _ 
dt  dxj 

(8) 

momentum  conservation: 

d{pUi)  ^  d(puiUj)  dp 

dt  '  dxj  dxi 

(9) 

energy  conservation: 

d{ph)  d{phuj)  dp  dp 

(10) 

dt  *  dxj  dt  ^  dxj 

Although  the  neglect  of  viscosity  is  a  great  simplification,  it  raises  interesting  problems 
for  the  modeling  of  those  flow  features  that  are  apparently  inviscid  but  require  some  viscosity 
for  initial  setup.  For  instance,  an  “exact”  solution  of  the  Euler  equations  requires  a  separate 
enforcement  of  the  Kutta  condition  in  order  to  avoid  infinite  pressures  at  the  trailing  edge. 
However,  many  difference  solutions  do  not  require  an  explicit  Kutta  condition  specification 
because  of  artificial  viscosity,  which  is  a  natural  result  of  discretization  errors.  A  hypothetical 
numerical  method  that  had  a  very  low  or  nonexistent  discretization  error  would  require  an 
explicit  enforcement  of  pressure  continuity  at  the  trailing  edge.  This,  however,  would  be  a 
very  small  price  to  pay  for  something  which  would  provide  great  benefits.  The  greatest  of 
these  benefits  would  be  the  ability  to  convect  vorticity  without  artificial  dissipation.  This 
ability  is  an  important  issue  for  rotorcraft,  on  which  the  blades  often  interact  with  vortices 
that  are  fairly  old,  but  have  not  dissipated.  By  contrast,  current  numerical  Euler  methods 
dissipate  vortices  in  a  very  short  time  (often  in  several  chords  of  travel). 

It  is  also  worth  noting  that  an  “exact”  inviscid  solution,  that  is,  one  lacking  dissipation, 
would  require  an  imposition  of  jump  conditions  in  order  to  predict  flows  with  shocks.  However, 
a  difference  solution,  with  proper  numerical  dissipation,  mimics  reality  and  permits  a  shock 
to  arise  naturally.  This  is  usually  referred  to  as  “shock  capturing.”  It  is  this  ability  that  was 
primarily  responsible  for  the  rapid  growth  of  CFD  in  the  1970’s.  The  presence  of  artificial 
viscosity  in  inviscid  methods  is  thus  both  a  curse  and  a  blessing. 

POTENTIAL  EQUATIONS 

A  considerable  simplification  of  the  inviscid  flow  equations  is  achieved  by  assuming  irro- 
tational  flow.  That  is, 

=  V  X  tt  =  0  (11) 

where  u  is  the  flow  vorticity.  Any  vector  field  «  that  is  irrotational  can  be  described  as  the 
gradient  of  a  scalar  potential  field  0,  or 


«  =  V<?i  (12) 

The  subsequent  replacement  of  three  velocity  components  by  one  scalar  potential  is  a  very 
useful  simplification.  However,  the  greatest  simplification  is  a  considerable  reduction  in  the 
number  of  equations  to  be  solved.  Consider  Crocco’s  equation, 

S  +  M  X  =  Vfio  d"  (13) 
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where  ho  is  the  stagnation  enthalpy.  (This  equation  is  derived  from  the  second  law  of  ther¬ 
modynamics  and  conservation  of  mass  and  momentum.)  If  a  flow  is  isentropic  (a  statement 
of  energy  conservation)  and  we  introduce  a  velocity  potential,  Eq.  (13)  reduces  to 


d(f> 

dt 

—  /iq  —  ) 

(14) 

Introducing  the  isentropic  relation 

II 

(15) 

we  obtain  Bernoulli’s  equation. 

pIpoo  =  1 

(16) 

which  together  with  the  continuity  equation, 

+  =  0  (17) 

forms  a  complete  equation  system.  The  greatest  simplification  of  the  potentizd  assumption  is 
that  it  permits  the  replacement  of  four  differential  equations  (conservation  of  momentum  and 
energy)  with  a  single  algebraic  equation.  Therefore,  potential  methods  are  a  minimum  of  five 
times  faster  than  equivalent  Euler  methods. 

Equations  (16)  and  (17)  are  actually  not  a  complete  system  for  lifting  solutions,  because 
the  presence  of  circulation  requires  a  branch  cut  so  that  the  potential  will  remain  a  single  valued 
function.  That  is,  it  is  necessary  to  specify  a  jump  in  potential,  F  =  A0  =  0u  —  0<,  whose 
magnitude  is  the  total  upstream  flow  circulation  (both  bound  and  free).  For  a  flow  with  free 
circulation  (vortex  and  circulation  sheets),  a  F  surface  must  coincide  with  the  vortex  sheets, 
which  usually  commence  at  trailing  edges.  An  expression  for  the  branch-cut  distribution  is 
obtained  by  enforcing  the  continuity  of  pressure  along  a  stream  surface  shed  from  the  trailing 
edge.  Using  Eqs.  (16)  and  (15),  we  have 


Ap  =  pu  -  p/  =  KA 


<t>t  +  -V0  •  V0 


=  0 


where  iF  is  a  constant,  or 


Ft  -|-  ave  ■  —  0 


(18) 


where  Kave  =  +  ^«)  is  the  average  velocity  of  the  circulation  sheet.  Equation  (18) 

expresses  the  Kutta  condition  and  it  completes  the  potential-flow  formulation.  This  Kutta 
equation  is  a  convection  expression  for  the  branch-cut  surface.  It  describes  the  surface  that 
shed  circulation  must  follow  in  order  that  there  are  no  forces  in  the  flow  field.  Equations  (18) 
and  (16)  are  thus  expressions  of  the  conservation  of  momentum.  Although  the  potential 
equat’  s  do  conserve  momentum  in  smooth  flows,  there  is  no  strictly  conservative  statement 
of  momentum  conservation.  The  resulting  jump  relations  permit  momentum  losses  through 
shocks,  and  constitute  a  nonphysical  drag  mechanism  in  the  isentropic  potential  model.  (The 
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physical  shock  drag  mechanism  is  a  jump  in  entropy  rather  than  a  momentum  across  a  shock.) 
This  physical  modeling  error  is  usually  small  and  easily  corrected.  Because  Eqs.  (16-18) 
strictly  conserve  mass,  they  are  referred  to  as  the  “conservative  full-potential  equations.”  With 
appropriate  treatment  of  shocks  and  vorticity  transport,  these  equations  probably  represent 
the  optimum  combination  of  accuracy  and  efficiency  for  inviscid  transonic  flow  problems. 


There  are  many  other  forms  of  the  potential  flow  equations.  Until  recently,  the  more 
usual  form  of  the  potential  equation  was  that  obtained  by  the  elimination  of  p  from  Eqs.  (16) 
and  (17).  That  is. 


6^0 


=  0 


(19) 


where  V  =  V(/>,  and  a  is  the  local  speed  of  sound  obtained  by  Bernoulli’s  equation  and  the 
isentropic  relation  =  Kp^~^.  Equation  (19)  is  readily  solved  and  forms  the  basis  of  severed 
widely  used  transonic  flow  codes.  Unfortunately,  Eq.  (19)  is  a  nonconservative  equation;  for 
transonic  methods  using  “shock  capturing,”  it  will  not  conserve  mass.  This  is  a  problem  of 
mathematical  consistency  for  which  there  is  no  simple  fix.  Furthermore,  there  is  no  significant 
cost  savings  in  using  the  nonconservative  mass  equation. 


The  above  potential  equations  (both  the  conservative  and  nonconservative  variants)  are 
usually  referred  to  as  “full-potential  equations”  because  no  terms  are  neglected  in  their  for¬ 
mulation.  However,  it  is  well  known  that  many  of  the  quadratic  and  higher  order  terms  in 
Eq.  (19)  must  be  very  small,  especially  for  flow  over  slender  bodies.  But  these  cannot  all  be 
neglected  for  it  is  among  these  terms  that  the  transonic  nonlinearity  is  to  be  found.  If  we 
neglect  all  the  higher-order  terms,  with  the  exception  of  the  quadratic  terms  in  the  streamwise 
direction,  we  obtain  the  simplest  approximate  equations  capable  of  predicting  transonic  flows. 
These  are  the  transonic  small-disturbance  equations,  of  which  there  are  many  variants,  all 
having  the  form 


A())tt  +  B(}>xt  =  Fx  +  4>zz  +  C(l)yy  -h  cross-  and  lower-order  derivatives  (20) 


where  the  nonlinearity  is  contained  in  the  streamwise  flux  F,  which  has  the  form  a0i  4-  /?0x- 
Equation  (20)  is  not  unique,  and  various  modifications  and  additions  have  been  made  to 
improve  its  stability  and  accuracy.  This  class  of  equations  is  often  referred  to  as  the  “TSD” 
(transonic  small  disturbance)  equations.  An  interesting  and  useful  feature  of  Eq.  (20)  is  that 
it  is  in  conservation  forru,  which  is  a  fortuitous  outcome  of  the  small-disturbance  limiting 
process.  The  next  level  of  simplification  is  to  eliminate  all  higher-order  terms  in  Eq.  (19), 
which  yields 

(21) 

where  Ooo  is  the  undisturbed  speed  of  sound.  This  is  the  linear  compressible  potential  equa¬ 
tion.  A  simple  Galilean  transformation  of  Eq.  (21)  yields  the  Prandtl-Glauert  equation.  A 
transformation  to  a  rotating  frame  does  not  produce  such  a  simple  result.  Equation  (21)  has 
been  used  to  predict  subcritical  compressible  rotor  flows.  Of  course,  this  equation  is  also  the 
basis  of  many  acoustic  rotor  analyses. 

The  final  simplification  is  to  assume  an  infinite  speed  of  sound  to  obtain  Laplace’s 
equation, 

V20  =  0  (22) 
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which  describes  incompressible  flow.  To  date,  most  rotorcraft  computational  development 
work  has  centered  on  solving  Eq.  (22).  The  reason  for  this  preponderance  of  work  is  the 
importance  of  accurately  predicting  the  wakes  of  helicopter  rotors.  Although  helicopter  rotors 
operate  deep  in  the  compressible  flow  regime,  the  evolution  of  their  shed  wakes  is  basically 
an  incompressible  problem.  Another  reason  for  the  importance  of  the  Laplace  equation  to 
rotorcraft  is  that  it  is  solvable  by  boundary-integral  methods  and  can  permit  very  rapid 
solutions.  This  efficiency  is  vital  for  problems  wherein  the  wake  solution  must  be  iterated 
with  solutions  for  the  blade  near-field  aerodynamics  and  motion. 

III.  SOME  IMPLEMENTATIONS  OF  THE  FLOW  EQUATIONS 

The  various  numerical  methods  and  codes  used  to  solve  the  above  systems  of  equations 
constitute  an  enormous  field  that  cannot  be  extensively  treated  in  the  scope  of  this  paper. 
Rather,  this  paper  will  be  restricted  to  those  methods  that  are  representative  of  present  work 
or  have  high  potential  for  engineering  application. 

As  previously  mentioned,  rotorcraft  CFD  must  follow  a  different  path  from  its  fixed-wing 
counterpart  because  of  its  unique  set  of  problems.  The  most  obvious  difference  between  the  two 
areas  is  the  greater  importance  of  potential  flow  methods  for  rotorcraft.  One  reason  for  this  is 
the  intrinsic  importance  of  unsteady  aerodynamics  to  rotorcraft.  (Even  in  fixed-wing  unsteady 
work,  such  as  flutter  prediction,  one  finds  a  predominance  of  potential  methods.)  Another 
reason  is  the  intrinsic  importance  of  wake  problems  to  rotorcraft.  The  viscous  problems 
encountered  on  rotorcraft  are  as  important  as  (and  more  challenging  than)  those  on  fixed- 
v/ing  aircraft.  Nevertheless,  unlike  in  the  fixed-wing  field,  the  nonviscous-dominated  problems 
of  rotorcraft  a:e  not  yet  solved,  and  these  require  much  of  our  present  effort. 

A  typical  CFD  review  would  be  mainly  concerned  with  differential  methods  of  solution 
(i.e.,  finite-difference,  finite- volume,  and  finite-element  methods  for  solution  of  the  partial- 
differential  equations).  However,  the  importance  of  the  wake  problem  to  rotorcraft  and  the 
difficulty  in  solving  it  with  difference  methods  is  such  that  boundary-integral  methods  must 
alsc)  be  considered.  Therefore,  the  prese;it  section  is  divided  into  a  discussion  of  both  differ- 
entieJ  and  integral  methods.  Differential  methods  will  be  discussed  in  terms  of  their  most 
common  use,  the  solution  of  local  blade  aerodynamics.  Integral  methods  will  be  discussed 
in  terms  of  their  most  common  application,  the  solution  of  complex  configuration  and  wake 
aerodynamics. 

BOUNDARY-INTEGRAL  METHODS 

The  prediction  of  wake  behavior  is  one  of  the  oldest  and  most  important  problems  in 
rotor  aerodynamics,  and  in  many  respects,  it  remains  unsolved.  Most  rotor  wake  problems 
are  incompressible.  These  problems  are  also  inviscid,  to  the  extent  that  wake  regions  are 
confined  to  thin  layers  that  are  convected  by  the  global  flow  but  are  otherwise  little  affected 
by  it,  somewhat  like  a  noninteracting  boundary  layer. 

Traditionally,  wakes  and  configurations  have  been  treated  by  means  of  boundary-integral 
methods.  The  basis  of  these  methods  is  the  ability  to  superpose  the  various  elemental  so¬ 
lutions  of  the  Laplace  equation — the  source,  th-e  doublet,  and  the  vortex.  The  actual  choice 
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of  functions  is  arbitrary  and  their  distribution  is  chosen  so  as  to  satisfy  surface  boundary 
conditions  (no  flow-through)  and  the  Kutta  condition. 

In  recent  years,  Morino  [1]  has  generalized  the  above  approaches  and  derived  an  integral 
flo  representation  that  is  valid  everywhere  (including  the  body  interior). 

The  basis  for  this  approach  is  the  second  Green’s  identity 

=  II  (.f  -  ^5  (23) 

where  G  is  the  Green’s  function,  satisfying 


V‘^G  =  6{x-x^) 


and  is  the  familiar  unit  source  potential 


G  =  - 


1 

4-jTr 


(24) 


(25) 


where  r  =  |x  —  x*|.  Green’s  identity  can  be  written  for  any  point  in  the  flow  held  (including 
the  body  interior)  as 


=  /£  - 11^  (26) 


where 


£'(a;»)  =  1,  outside  Sb 

=  ^,  on  5b  (27) 

=  0  inside  Sb 


Note  that  the  body  is  represented  by  a  source  and  doublet  distribution  on  its  surface,  and 
the  wakeis  represented  by  a  doublet  layer.  If  x*  is  in  the  held  (outside  the  body),  Eq.  (27)  is 
an  integral  representation  for  0  in  terms  of  the  values  of  0  and  d(}>/dn  on  the  surface  of  the 
body,  5b,  and  of  F  on  the  surface  of  the  wake,  Sw-  If  is  on  the  surface  5b,  Eq.  (27)  is 
a  compatibility  condition  between  the  values  of  0  and  50/5n  on  5b,  and  of  F  on  Sw-  Note 
that  d(j)/dn  on  5b  is  kno./n  from  the  surface  boundary  condition,  and  that  F  on  Sw  is  known 
from  the  preceding  time  history.  Therefore,  Eq.  (27)  is  an  integral  equation  that  may  be  used 
to  directly  evaluate  0  on  5b. 

Once  0  on  the  body  surface  is  known,  Eq.  (27)  can  be  used  to  calculate  0,  and  hence  to 
calculate  the  velocity  and  pressure  anywhere  in  the  field.  In  particular. 


V-V*0 


=  // 

JJs,  [an 


:v.g-4(v.g) 


dS  -  J J r-^{V,G)dS  (28) 


This  equation  permits  us  to  calculate  the  velocity  of  the  wake  points  and  hence  the  geometry 
of  the  wake  at  a  time  t.  Note  from  Eq.  (18)  that  F  follows  these  wake  points.  Also,  since  the 
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Kutta  condition  of  pressure  continuity  implies  that  the  trailing-edge  circulation  is  shed  into 
this  wake  at  every  time  step,  we  have  conditions  that  determine  this  quantity  at  the  upstream 
edge  of  the  wake. 

The  above  equations  are  the  basis  for  all  boundary-integral  methods.  A  variant  of  this 
formulation  [2]  results  from  expressing  Eq.  (27)  at  the  inside  and  outside  surfaces  of  the  body 
and  combining  these  to  obtain 


where  Sb  ~  P  indicates  the  exclusion  of  point  P  in  the  first  integral.  Equation  29  gives  the 
total  potential  at  an  interior  point  i  as  the  sum  of  perturbation  potentials  that  result  from 
normal  doublet  distributions  of  strength,  {<f>o  —  4>i)  on  Sb  and  (0^  —  0/)  on  Sw,  and  from  a 
source  distribution  of  strength,  {d4>/dn)o  —  {d(i)jdn)i  on  Sb-  In  principle,  an  infinite  number 
of  combinations  of  doublet  and  source  distributions  will  give  the  same  external  flow  field,  but 
different  internal  flow  fields.  To  render  a  unique  combination  of  singularities,  we  can  specify 
one  of  the  singularity  distributions,  or  the  internal  flow.  Two  obvious  choices  of  internal  flows 
are  stagnation  and  undisturbed  flow. 

The  above  formulations  are  solved  by  discretizing  the  body  surface  and  wake  into  suitably 
conforming  elements  (or  panels)  and  expressing  the  above  equations  as  discrete  summations. 
The  resulting  system  of  linear  equations  constitutes  a  large,  full  matrix  which  is  readily  solved 
by  direct  or  iterative  methods.  These  methods  are  usually  referred  to  as  “panel  methods.” 
Panel  methods  vary  according  to  the  type  of  elemental  solutions  used,  the  manner  of  dis¬ 
tributing  them,  and  the  type  of  discretization.  The  most  common  methods  in  use  are  “low 
order”  schemes,  in  which  the  singularity  distribution  is  assumed  to  be  constant  on  each  panel. 
“High  order”  methods  assume  a  variation  of  source  or  doublet  strength  over  each  panel. 
Some  schemes  entirely  ignore  thickness  effects  and  represent  wings  or  rotors  as  singularity 
planes  with  no  source  terms.  These  methods  are  referred  to  as  “lifting  surface  methods.”  An 
even  greater  simplification  involves  the  approximation  of  this  surface  as  a  single  streamwise 
element — that  is,  by  “lifting  line  methods.”  The  above  equations  represent  solutions  in  terms 
of  source  and  doublet  distributions.  A  lattice  of  line  vortices  can  be  shown  to  be  equivalent 
to  an  array  of  constant-strength  doublet  panels  defined  on  the  same  vortex  lattice.  In  “vor¬ 
tex  lattice  methods,”  velocities  are  computed  using  the  Biot-Savart  law.  The  majority  of 
rotor/wake  analyses  are  lifting-line  and  vortex-lattice  methods. 

Most  recent  boundary-integral  methods  have  involved  the  use  of  panel  methods.  The 
greatest  impetus  for  these  developments  has  been  the  need  for  modeling  complex  aircraft 
(usually  fixed-wing)  geometries.  Perhaps  the  most  widely  used  panel  method  codes,  called 
“VS AERO”  [3,  4]  and  “PAN AIR”  [5,  6],  are  based  on  source/doublet  representations  such  as 
Eq.  (29).  VS  AERO,  in  particular,  has  found  wide  use  for  rotorcraft  component  interaction 
problems,  because  it  models  .separation  effects  and  free-wake  convection.  Figure  2  shows 
typical  VSAERO  computations  of  a  complex  rotorcraft  configuration  [7].  The  free  convection 
of  the  wing  and  rotor  wakes  of  a  V-22  can  be  clearly  seen.  A  comparison  of  PANAIR  and 
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VSAERO  for  a  V-22  in  airplane  mode  is  made  in  Ref.  [8].  This  latest  generation  of  integral 
codes  has  been  driven  mainly  by  fixed-wing  needs.  Rotor  hover  problems,  however,  have  not 
been  extensively  treated  by  panel  methods. 

Most  panel  method  codes  have  been  designed  for  fixed-wing  applications  and  require 
only  a  specified  wake  geometry.  However,  for  low-speed  flight,  the  wake  induction  becomes 
increasingly  important  compared  to  free-stream  convection.  For  hovering  helicopters  this 
wake  induction  is  the  dominant  flow  consideration,  especially  because  the  rotor  remains  close 
to  its  own  wake.  In  fact,  helicopter  wake  problems  have  mainly  been  treated  by  vortex  lattice 
methods  using  a  lifting-line  representation  of  the  blade. 

The  hover  problem  is  very  sensitive  to  the  accuracy  of  the  wake  computation.  This  places 
a  stringent  requirement  on  the  wake  panel  or  element  modeling  density.  In  order  to  reduce  the 
size  of  the  wake  grid,  various  efforts  have  been  made  to  use  higher-order  geometry  definitions 
and/or  integration  of  the  wake  elements.  The  most  current  example  of  this  is  the  curved 
vortex  elements  of  Bliss  [9,  10].  In  this  work  every  three  wake  points  on  a  vortex  line  are 
used  to  define  a  locally  parabolic-shaped  element  whose  induced  velocity  is  described  by  an 
analytic  expression  (Fig.  3).  The  accuracy  of  this  element  is  such  that  far  fewer  points  are 
required  to  resolve  the  wake  than  if  a  series  of  straight-line  elements  were  used. 

The  physical  accuracy  of  the  hover  wake  model  is  also  a  critical  item.  One  physical 
wake  feature  that  is  uniquely  important  to  rotor  hover  wake  modeling  is  the  local  self-induced 
velocity  that  results  from  wake  curvature.  A  curved-line  vortex  segment  contains  a  logarithmic 
singularity  that  does  not  occur  in  reality  because  real  vortices  have  non-zero  core  radii.  But  the 
wake  settling  rates  are  a  function  of  the  core  size.  Although  this  is  a  logarithmic  function  (i.e., 
not  strong)  it  is  important  for  hover.  It  has  been  shown  [11,  12]  that  the  effect  of  this  core  size 
can  be  accounted  for  by  incorporating  an  appropriate  cutoff  distance  in  the  local  Biot-Savart 
integration  over  the  curved  element  (Fig.  4).  Most  wake  analyses  use  straight-line  vortex 
elements  only.  The  contribution  at  any  wake  collocation  point  from  its  two  adjacent  vortex 
elements  is  therefore  zero.  This  corresponds  to  a  cutoff  distance  of  one  wake-element  interval, 
and  is  too  large.  Such  refinements  are  not  required  for  computations  for  many  advancing 
rotors  or  fixed-wings,  because  the  free-stream  convection  overwhelms  these  effects. 

Although  free- wake  computations  for  hovering  rotors  have  been  performed  for  many  years, 
the  most  reliable  computations  have  tended  to  be  those  that  use  empirical  wake-geometry 
data  [13,  14,  15].  These  codes  are  also  extremely  fast,  because  they  do  not  require  inflow 
computations  on  the  wake.  Such  codes  can  produce  good  thrust/power  polars  when  a  suitable 
wake  data  base  is  available  (a  good  airfoil  data  base  is  also  required).  However,  there  is 
a  fortuitous  element  in  this  approach,  since  these  codes  often  do  not  predict  the  thrust- 
vs.-collective-pitch  curves  accurately.  The  earlier  free-wake  codes  did  not  predict  the  wake 
geometry  well  and  did  not  give  good  performance  predictions.  More  recent  free-wake  codes 
are  now  beginning  to  produce  good  results,  because  of  better  numerical  techniques  that  more 
accurately  model  the  physics.  The  code  EHPIC  [16]  is  one  that  is  producing  good  wake  and 
performance  predictions.  It  is  a  vortex-lattice,  lifting-surface  code  that  uses  the  self-induction 
model,  curved  wake  elements,  and  a  unique  wake-relaxation  scheme,  which  obviates  many 
instabilities. 

A  recent  paper  by  Felker  et  al.  [17]  gives  good  examples  of  current  free-wake  capabilities 
using  the  EHPIC  code  and  data  from  many  tested  rotors.  Figure  5  shows  a  fairly  typical 
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comparison  of  measured  and  computed  tip- vortex  geometries  and  lift/power  polars  for  an 
S-76  rotor.  The  comparison  of  tip-vortex  axial  settling  rate  shows  some  differences, but  the 
distance  of  closest  blade/vortex  approach  is  close.  Furthermore,  the  predicted  and  measured 
power  are  within  3%  of  each  other.  This  example  is  shown  here  only  because  it  involves  both 
wake  and  performance  data  for  a  modern  rotor.  The  level  of  agreement  between  data  and 
computation  is  good.  Another  widely  used  code  is  HOVER  [18],  which  also  uses  a  lifting- 
surface,  free- wake  analysis.  HOVER  also  employs  a  local  wake  curvature  treatment. 

Vortex-lattice  codes  are  also  widely  used  for  forward-flight  wake  analyses.  The  Scully 
vortex-lattice  model  [19]  is  incorporated  in  the  CAMRAD  comprehensive  rotor  analysis  code 
[20]  and  is  probably  the  most  widely  used  of  such  codes.  The  model  uses  a  network  of 
streamwise  vortices  (generated  by  the  bound-circulation,  spanwise  gradient)  and  spanwise 
vortices  (generated  by  the  time  rate  of  change  of  bound  circulation).  An  alternate  vortex- 
lattice  model  has  also  been  developed  that  uses  the  Bliss-curved  vortex  elements  (Fig.  6). 

These  advancing  rotor  applications  tend  to  be  used  with  lifting-line  blade  models.  These 
blade  models  use  measured  airfoil  data  to  obtain  profile  drag  and  pitching  moments.  However, 
this  approach  precludes  the  treatment  of  3-D,  unsteady  transonic  characteristics.  It  also 
requires  empirical  models  for  stall.  More  accurate  treatment  of  the  blades  requires  the  use  of 
differential  CFD  methods. 

DIFFERENCE  METHODS 

The  computation  of  the  rotor-blade  aerodynamics  entails  the  treatment  of  various  nonlin¬ 
ear  effects,  the  most  common  of  which  is  transonic  flow.  The  simplest  equations  for  treating 
transonic  flows  are  the  nonlinear  potential  equations.  Because  these  equations  are  nonlin¬ 
ear,  they  cannot  be  superposed  and  it  is  necessary  to  discretize  and  solve  them  directly.  We 
shall  discuss  this  process  only  for  potential  methods.  However,  the  basic  process  has  much  in 
common  with  that  for  the  more  complex  flow  approximations,  such  as  Euler  and  Reynolds- 
averaged  Navier-Stokes. 

The  first  step  for  helicopter  computations  is  to  transform  the  equations  to  a  translating 
and  rotating  coordinate  system  fixed  to  the  blade, 

=  Uoot  +  Dt  X  r  (30) 

where  r'  =  {x',  y',  z')  and  r  =  (x,  y,  z)  are  the  inertial  and  body-fixed  coordinates  (Fig.  7),  and 
UoQ  and  D  are  the  translational  and  angular  velocities  of  the  rotor.  Under  this  transformation, 
the  conservative  potential  equations,  Eqs.  (16)  and  (17),  become 

Pt+V-  \pV]  =  0  (31) 

and 

p/Poo  =  |l  -  ^  A  +  -  v^)  I’"'  (32) 

where 

V  =  Voo  +  V0  (33) 
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is  the  local  velocity  in  the  blade-fixed  frame  and  V^o  =  C^oo  +  ft  x  r  is  the  local  free-stream 
velocity  seen  by  an  observer  in  the  rotating  frame.  Of  course  0,  being  scalar,  is  imchanged 
by  the  coordinate  transformation.  But  in  the  blade-fixed  frame  (Eq.  (33))  it  appears  as 
a  perturbation  about  the  free-stream  velocity,  Voo-  Note  that,  since  Vco  is  rotational,  we 
cannot  combine  it  into  (or  define)  a  “full  potential”  in  the  normal  fixed-wing  sense. 

Equations  (16)  and  (17)  are  commonly  written  for  a  generalized,  moving  coordinate  sys¬ 
tem  T],  C,  t  (fig.  8)  as 

dtip/J)  +  d^ipU/J)  +  dr,{pV/J)  +  d^ipW/J)  =  0  (34) 

where  U,  V,  and  W  are  contravariant  velocity  components,  and  J  is  the  Jacobian  of  the 
coordinate  transformation.  Bernoulli’s  equation  is  similarly  expressed  in  general  coordinates: 

p/poc  =  {l  -  ^  [4>t  +  {U  +  +  iV  +  vt)cl>r,  +  iW  +  O)0c]} ""  (35) 

L  ^oo  J 

In  this  formulation,  the  blade  motion  is  specified  by  the  coordinate  terms  Ct-  Equations 
(34)  and  (35)  represent  an  excellent  formulation  of  the  potential  problem  because  they  are 
conservative,  compact,  and  geometrically  general. 

By  contrast,  under  a  rotational  and  translational  transformation,  the  nonconservative  full- 
potential  equation  (Eq.  (19))  becomes  (in  Cartesian  coordinates  for  a  rotor  in  pure  edgewise 
motion) 


<j)tt  4-  2Vx4>xt  +  2Vz4>zt  =  —  V^)<i>xx  +  —  Vz)4>zz  +  (o^  ~ 

-  2VxVy(t>xy  -  2VxVz(f>xz  -  2VyVz(l)yz  (36) 

-I-  (ft^x  —  2QUoo  cos  ^)0i  +  (ft^^  -I-  2ftt/oo  sin  i>)(j>z 

which  cannot  be  put  in  conservation  form.  The  small-disturbance  equation  retains  the  form 
shown  in  Eq.  (20)  under  a  rotational  transformation. 

The  first  differential  CFD  methods  developed  were  mainly  concerned  with  the  small- 
disturbance  equation  and  much  rotor  work  is  still  done  with  these  techniques.  This  work  uses 
finite  difference  methods,  in  which  the  partial  derivatives  in  the  basic  equations  are  directly 
approximated  by  differences  in  order  to  derive  systems  of  linear  equations.  Other  discretization 
methods  have  since  been  developed,  notably  finite- volume  and  finite-element  methods.  These 
methods  also  work  well,  the  finite-element  methods  being  the  most  recent.  The  following 
discussion  will  focus  on  the  finite-difference  approach. 

In  order  to  illustrate  basic  ideas  and  issues,  we  will  consider  a  simple  problem  which 
is  a  paradigm  for  all  finite-difference,  potential  methods.  Consider  a  simple  2-D  transonic 
small-disturbance  equation  (see  Fig.  9) 


in  which 


/50XI  d"  0yy  —  0it 


0  =  + 


(37) 
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This  equation  models  low-frequency  flows  that  result  in  large  shock  hysteresis.  To  complete 
this  problem  we  require  boundary  conditions  and  the  Kutta  condition  for  pressure  continuity. 
The  latter  is  satisfied  by  convecting  the  trailing-edge  circulation  using 

Tt  +  =  0  (38) 


The  body  no-flow-through  condition  is  specified  on  a  mean  surface  (y  =  0)  and  therefore  takes 
a  transpiration  form 

*^2/u  ~^oo/u(^i^)  (39a) 

on  the  upper  side  of  the  mean  surface  and 


^yt  —  ^<x>fi  (^)  0 


(396) 


on  the  lower  side  of  the  mean  surface,  where  the  body  geometry  is  given  by  y  =  /u(x),/!(x) 
for  the  upper  and  lower  surfaces.  Various  far-field  boundary  conditions  can  be  used.  For  an 
unsteady  problem  the  best  condition  to  use  is  a  nonreflecting,  first-order  wave  equation  which 
permits  no  inward  propagation  from  the  far-field  boundary.  However,  the  specification  of 
nonperturbed  flow  is  very  common  and  simple.  Therefore  we  will  specify  that  at  the  far-field 
boundary 

(f>  =  UooX  (40) 

In  order  to  solve  this  system  of  equations  (Eqs.  (37-40))  we  construct  a  Cartesian  grid  and 
approximate  the  various  partial  derivatives  by  differences.  The  following  difference  approxi¬ 
mations  will  now  be  applied  to  Eq.  (37): 


1st  backward  difference  in  x,  =  Vx(f>  =  {<i>i  —  4>i-\)IDx 
1st  forward  difference  in  x,  <t>x  =  Ay0  =  —  (l>i)IDx 

2nd  backward  difference  in  x,  (i>xx  —  VxVx0  =  (0i  —  20i_i  -f  (l)i-2)/Dx^ 
2nd  centered  difference  in  x,  <l>xx  —  ViAx0  =  (0i-i  —  20,  -I-  0i_i)/£)x^ 


(41) 


where  Dx  is  a  streamwise  grid  interval.  These  expressions  are  easily  derived  from  Taylor 
series.  Modifications  of  these  difference  approximations  are  required  for  points  adjacent  to  the 
boundaries.  For  instance,  for  points  adjacent  to  the  mean  airfoil  surface  the  normal  second 
derivatives  are  expressed  as 


^yyj  ~  i^y^J  ^yn)l^y 

=  (0y+i  -  0J  -  DyUoofi{x))IDy‘^ 

(42a) 

on  the  upper  surface,  and 

<t>yyj  =  {DyUooflix)  -  0j  +  (f>j-i)/Dy^ 

(426) 

on  the  lower  surface.  The  bound  and  shed  circulation  of  a  lifting  problem  require  branch  cuts 
(potential  jumps  or  discontinuities)  in  the  wake.  Adjacent  to  these  branch  cuts,  the  normal 
second  derivatives  are  expressed  as 

^yyj  =  {^y4>J  “  '^y<t>j)IDy 

=  -  2<t>J  +  0J-1  -  ^)/Dy'^ 

(43a) 
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above  the  cut,  and 


(t>yyj  =  -  r)/Dy2  (436) 

below  the  cut.  These  expressions  for  differencing  across  the  cut  are  easily  obtained  by  com¬ 
bining  a  Taylor  expansion  for  <p  with  the  known  discontinuities  in  <j>  and  d4>jdn  (F  and  0, 
respectively).  The  value  of  Fj  is  found  by  solving  a  discretized  form  of  Eq.  (38). 

An  instructive  (but  inadvisable,  for  stability  reasons)  discretized  form  of  Eq.  (38)  is 


rr^‘  -  rr  .  „  t?  -  r?-,  „ 

Af  Ax  “ 

which  is  solved  by  a  simple  downstream  marching  process, 

r;‘+‘  =  (i-a)r?+Qr;'_i 


(44) 

(45) 


where  a  =  C/ooAt/Aa:.  Equation  (44)  evaluates  the  spacial  differences  at  the  previous  time 
step,  n,  and  F"+^  is  solved  as  a  simple  algebraic  function  of  known  quantities  at  the  previous 
time  step.  Such  a  scheme  is  referred  to  as  an  “explicit”  method.  For  a  >  1,  it  is  unstable. 
Note  that  for  a  =  1  we  have  F”'*'^  =  Ff_i.  This  particular  choice  of  grid  amounts  to  physically 
convecting  individual  fluid  particles  from  point  to  point,  which  is  Lagrangian  convection,  and 
numerical  dissipation  does  not  occur.  In  general,  however,  grids  cannot  be  so  chosen  and 
dissipation  does  occur.  Such  dissipation  is  unimportant  for  problems  in  which  circulation 
merely  convects  away  from  the  body,  as  in  this  example.  However,  this  is  an  important 
matter  for  recirculant  flows  such  as  the  hover  problems  mentioned  in  the  previous  section. 
There,  the  rotor  boundary-integral  methods  treated  the  wake  by  Lagrangian  convection.  This 
type  of  convection  is  not  done  with  current  differential  CFD  methods  (potential,  Euler,  etc.). 
An  additional  feature  of  potential- flow  CFD  is  that  the  wake  differencing  (such  as  in  Eq.  (43), 
which  is  the  universal  approach)  confines  the  wake  to  a  grid  plane.  (We  shall  later  show  that 
this  is  not  necessary.)  By  contrast,  a  conventional  Euler  solver,  having  no  potential  jump, 
convects  vorticity  freely  and  naturally;  but  it  still  does  not  solve  the  dissipation  problem. 

In  our  simple  problem,  we  first  consider  a  subsonic  flow — that  is,  >  0.  A  discretization 
of  Eq.  (37)  that  is  suitable  for  subsonic  flow  is 

-  h  +  Vj,Aj,0"+i]  =  (46) 

where  h  =  Di.  Note  that  the  entire  array  of  unknowns  must  be  found  simultaneously  in 
a  large  matrix  inversion.  Such  a  scheme  is  called  “implicit”  and  has  the  advantage  of  having 
no  linear  stability  restrictions  (unlike  the  above  explicit  method).  However,  the  price  for  this 
stability  is  the  inversion.  The  left-hand  side  is  nearly  Laplacian,  and  can  be  inverted  iteratively 
by  various  standard  methods.  Furthermore,  since  and  0”  are  close,  the  iteration  process 
must  be  very  fast  (requiring  as  little  as  one  iteration  for  some  methods).  However,  a  faster 
and  more  usual  approach  involves  an  approximate  factorization  of  the  left-hand  operator  into 
a  more  easily  inverted  form.  For  Eq.  (46)  such  a  factorization  is 

(/  -  /3hAx)(Vx  -  /iV„Ai,)0"+^  =  Vx0”  -I-  /?h2A:,VyAy0"  (47) 
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The  new  right-hand  term  is  the  factorization  error,  which  is  evaluated  explicitly  at  the 
previous  time  step.  Equation  (47)  is  very  easily  solved  as  a  two-step  process: 

(7  -  =  Vx</>"  +  /3h^A^VyAy(l>^  (48a) 

(Vx  -  hVyAy)<j)^+^  =  0*  (486) 

First  the  upper  bidiagonal  inversions  (Eq.  (48a))  are  performed  on  each  y=constant  grid 
line.  Then  the  lower  triangular  matrix  (Eq.  (486))  is  solved  on  each  x=constant  grid  line. 
These  are  very  simple  and  well-behaved  processes.  This  often-used  procedure  is  referred  to  as 
“approximate  factorization”  (AF).  Many  types  of  AF  are  possible,  and  it  can  be  shown  that 
the  well-known  ADI  (alternating  direction  implicit)  method  is  a  form  of  AF.  For  supersonic 
flows  (when  /?  <  0)  the  difference  equation  (Eq.  (46))  is  unstable,  and  it  becomes  necessary 
to  replace  the  streamwise  term  by  a  backward  difference,  such  as  <^xx  —  VxVx0.  The  first 
successful  transonic  potential  codes  employed  a  simple  switch  in  the  differencing  process  (based 
on  the  sign  of  0)  to  treat  mixed  supersonic-subsonic  flows.  Almost  all  subsequent  transonic 
algorithms  have  employed  some  related  forms  of  upwind  biasing. 

This  discussion  is  meant  only  to  expose  the  main  issues  and  must  leave  out  many  im¬ 
portant  related  issues.  For  instance,  the  above  discussion  uses  a  nonconservative  equation  for 
simplicity.  However,  the  basic  equations  and  their  difference  forms  must  remain  in  conser¬ 
vation  form  in  order  to  properly  capture  shocks.  The  construction  of  good  algorithms  is  a 
very  extensive  field,  requiring  much  work  on  differencing,  stability,  accuracy  (including  con¬ 
servation),  inversion  methods,  boundary  conditions  and  grid  generation.  A  review  of  these  is 
beyond  the  scope  of  this  paper.  However,  this  example  has  exposed  some  of  the  salient  issues 
and  terminology  for  these  solution  methods.  With  this  background,  we  will  now  discuss  some 
of  the  actual  codes  that  have  been  developed. 

Because  Eqs.  (20)  (a  TSD  equation)  and  (36)  (a  nonconservative  full-potential  equation) 
have  the  same  general  form,  they  share  the  same  general  solution  methods.  Equation  (20)  is 
stably  discretized  using  a  “Murman  mixed  difference  scheme”  (i.e.,  the  previously  mentioned 
central/backward  differences  in  the  subsonic/supersonic  flow  regions  [21])  in  the  streamwise 
direction.  The  counterpart  to  the  Murman  scheme  for  Eq.  (36)  is  the  “Jameson  rotated 
method,”  in  which  mixed  differences  in  all  directions  are  used  [21].  Equation  (20)  has  been 
implemented  in  the  finite-difference  rotor  code  (FDR)  [22,  23]  (which,  in  spite  of  the  name, 
is  conservative)  which  uses  an  ADI  solution  scheme.  A  refined  version  of  this  code  called 
TSP  [24]  is  the  best-developed  and  most  heavily  used  small-perturbation  code  for  rotors. 
The  nonconservative,  full-potential  equation  (Eq.  (36))  has  been  implemented  in  several  rotor 
codes.  The  first  such  implementation  is  the  steady  code  ROT22  [25] — a  derivative  of  FL022, 
which  uses  a  successive-line  over-relaxation  (SLOR)  inversion  scheme.  Another  steady  code 
for  Eq.  (36)  is  TFARl  [26],  which  uses  an  approximate  factorization  solution  method  (rather 
than  SLOR).  The  fully  unsteady  form  of  Eq.  (36)  is  implemented  in  the  code  TFAR2  [27],  a 
derivative  of  TFARl. 

The  above  codes  are  limited  either  by  small-disturbance  or  conservation  considerations. 
The  implementation  of  a  conservative  full-potential  method  is  complicated  by  the  inability  to 
either  directly  eliminate  p  (which  destroys  the  conservation  form)  or  to  stably  solve  Eqs.  (34) 
and  (35)  as  a  two-equation  system.  The  spacial  terms  in  Eq.  (34)  present  no  problem,  as  these 
are  easily  treated  using  central  differences  (stability  being  maintained  by  upstream  biasing  of 
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density,  with  density  evaluated  at  the  previous  time  step  or  iteration).  The  real  problem  is  to 
express  the  term  pt  in  a  form  that  does  not  depend  on  p  at  the  new  time  step,  n  +  1.  This 
can  be  achieved  [28,  29]  by  expanding  p^'^^  as 

(r+'  -  4,--)  (49) 

where  dpjd(t>  is  a  differential  operator  obtained  from  Bernoulli’s  equation.  When  Eq.  (49)  is 
substituted  into  the  discretized  form  of  Eq.  (34)  there  results  an  equation  of  the  form 

+  V5(p”A44>"+‘)  +  +  Vc(4."Acr+‘)  =  C" 

This  equation  can  be  considered  to  be  a  hybrid  of  Eqs.  (34)  and  (36).  The  spacial  flux 
operators  correspond  to  those  in  the  conservative  full-potential  equation,  Eq.  (34).  However, 
the  time  derivatives  have  the  same  form  as  the  nonconservative  time  terms  of  Eq.  (36)  do.  The 
resulting  conservation  error  is  cancelled  by  the  term  C”,  which  is  a  conservation  correction 
term  evaluated  at  time  n.  This  formulation  was  first  implemented  in  a  fixed-wing  code  called 
TUNA  [29].  A  subsequent  derivative  of  this  code  called  FPR  [30,  31]  is  tailored  for  rotor 
applications.  These  codes  both  use  an  approximate  factorization  inversion  scheme.  A  similar 
rotor  implementation  called  RFS2  [32]  uses  a  strongly  implicit  procedure  for  matrix  inversion. 

All  of  the  above-mentioned  rotor  codes  (ROT22,  FDR,  TSD,  TFARl,  TFAR2,  FPR, 
RFS2)  have  been  used  in  industry.  ROT22  is  among  the  most-used  because  its  simplicity 
and  robustness  do  not  put  great  demands  on  the  user.  This  code  is  probably  best  used 
for  initial  high-speed-configuration  comparisons  (those  based  on  planform  and  profile,  for 
instance).  However,  its  shock  errors  (a  result  of  being  nonconservative)  and  its  inabifity  to 
handle  unsteady  effects  limit  its  accuracy  and  load-prediction  ability.  In  the  U.S.,  FPR  is 
probably  the  most  available  and  highly  developed,  both  from  a  technical  and  from  a  user 
viewpoint,  and  it  has  a  sizable  user  community. 

At  this  point  it  is  useful  to  demonstrate  some  of  the  capabilities  of  differential  CFD 
methods  to  predict  the  local  blade  flows  on  helicopter  rotors.  This  will  be  done  by  comparison 
of  computed  flows  with  results  from  experiments  that  were  designed  for  this  purpose. 

Surface  pressure  data  from  nonlifting  rotors  validate  the  basic  ability  of  codes  to  predict 
transonic  flows,  free  from  the  complications  of  wake  effects.  These  tests  are  usually  performed 
on  advancing  rotors  in  order  to  avoid  the  wake  buildup  that  would  occur  on  a  nonlifting  hover. 

Probably  the  most  extensive  surface  pressure  model  rotor  data  base  has  been  acquired  by 
ONERA  [33,  34,  35].  This  data  includes  nonlifting  and  lifting  data  for  two-  and  three-bladed 
rotors  with  a  variety  of  blade  profiles  and  planforms.  Figure  10  shows  an  early  nonUfting 
computation  of  the  pressure  variation  on  two  surface-pressure  transducers  performed  with  a 
2-D  small-disturbance  code.  Subsequent  3-D  computations  compare  equally  well  with  the 
data.  However,  the  point  of  discussing  an  inboard  2-D  computation  is  to  show  the  importance 
of  transonic  unsteadiness  without  any  m'tigating  3-D  influences.  This  unsteadiness  is  seen 
in  the  asymmetry  of  the  pressure  about  the  rp  =  90°  azimuth.  Computations  for  a  higher 
speed  case  in  which  the  steady  and  unsteady  full-potential  codes  TFARl  and  TFAR2  [27] 
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were  used  are  shown  in  Fig.  11.  The  inadequacy  of  steady  computations  is  clearly  shown  in 
this  comparison. 

Figures  12  and  13  show  a  nonlifting  computation  and  data  comparison  [36]  that  is  unique 
in  that  it  employs  optical  methods  rather  than  the  usual  pressure  instrumentation.  Figure 
12  shows  a  series  of  inter ferograms  which  are  obtained  from  holograms  produced  for  a  range 
of  rotor/laser-beam  orientations.  These  interferograms  can  be  used  to  reconstruct  nearly  the 
entire  pressure  field  around  the  rotor.  Figure  13  shows  these  reconstructed  pressure  fields  at 
0.08  chords  above  the  rotor  and  at  several  radii.  Also  included  in  Fig.  13  are  a  comparison 
between  small-disturbance  computations  and  laser  velocimetry  results.  Although  several  un¬ 
explained  differences  are  seen,  the  overall  comparison  is  quite  good.  Such  comparisons  have 
demonstrated  the  ability  of  CFD  methods  to  predict  pressure  fields  away  from  the  surface  of 
the  rotor — an  essential  for  high-speed  acoustics. 

Nonlifting  experimental  data  have  also  been  used  to  validate  the  ability  of  codes  to  predict 
high-speed-profile  power  and  transonic-drag  rise.  Figure  14  shows  a  comparison  of  computed 
and  measured  torque  for  a  nonlifting  two-bladed  rotor  with  a  NACA0012  airfoil.  The  computa¬ 
tion  [37]  was  performed  using  both  the  standard  FPR  code  and  a  variant  with  a  nonisentropic 
correction.  This  computation  also  used  a  Nash-MacDonaJd  boundary  layer  model.  The  torque 
was  obtained  by  integrating  the  skin-friction  and  surface  pressures.  The  transonic  drag  diver¬ 
gence  is  well  predicted,  whereas  the  purely  isentropic  computation  overpredicts  the  torque  at 
the  highest  tip  Mach  Numbers. 

These  and  many  other  related  validations  have  demonstrated  that  potential  methods  can 
predict  those  basic  blade  flows  that  are  essentially  inviscid.  This  means  they  can  predict 
nearly  all  flows  except  those  that  involve  retreating-blade  stall.  Much  remains  to  be  done 
on  these  problems,  especially  on  applying  boundary- layer  corrections  to  improve  drag-  and 
pitching- moment  calculations.  Another  concern  is  the  manner  of  combining  these  near-blade 
computations  with  the  wake-aerodynamic  models  to  obtain  a  global  rotor-flow  analysis. 

IV.  THE  COMPUTATION  OF  COMPLETE  ROTOR  FLOWS 

The  previous  section  summarized  some  of  the  last  twenty  years  of  CFD  code  development 
related  to  rotorcraft.  There  now  exists  an  extensive  arsenal  of  hard- won  computational  tools 
with  which  to  attack  rotor  flow  problems.  However,  rotor  behavior  involves  so  many  interre¬ 
lated  phenomena  that  the  effectiveness  of  any  one  analysis  can  be  difficult  to  judge.  The  need 
for  complete  analyses  is  obvious  and  has  been  addressed  elsewhere  in  the  context  of  comprehen¬ 
sive  code  development.  These  comprehensive  codes  have  used  a  variety  of  boundary-integral, 
analytical,  and  empirical  aerodynamic  methods.  An  important  requirement  for  such  methods 
is  that  they  be  fast  enough  to  be  an  integral  part  of  a  total  vehicle  analysis.  The  following 
discussion  addresses  the  capability  of  CFD  for  modeling  the  total  rotor  aerodynamics,  and  the 
methods  for  at  least  including  CFD  in  a  comprehensive  rotor  analysis.  The  first  part  of  this 
discussion  will  involve  the  development  of  hybrid  methods  in  which  various  wake  and  blade 
elements  are  assembled  to  produce  a  total  aerodynamic  analysis.  An  additional  section  will 
discuss  some  methods  for  a  unified  CFD  analysis  of  blade-wake  flow  systems. 
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HYBRID  METHODS 

Figure  15  illustrates  the  basic  computational  problem  for  a  hovering  rotor.  The  behavior 
of  a  hovering  rotor  is  governed  by  its  system  of  shed  vortices  and  vortex  sheets  and  by  the 
local  blade  flow.  Superimposed  on  the  blade  (Fig.  15)  is  a  representation  of  a  typical  grid 
used  for  computing  local  rotor  flows.  The  grid  extends  about  3-10  chordlengths  from  the 
blade  in  all  directions.  The  computation  of  the  flow  in  this  grid  differs  from  that  in  the  simple 
example  in  Section  III  in  that  many  of  the  shed  vortex  and  sheet  elements  pass  through  this 
grid  and  come  physically  close  to  the  blade.  These  circulatory  elements  must  be  inserted  into 
the  near- held  computation.  After  this  is  done,  we  still  have  an  incomplete  problem  with  this 
type  of  blade-oriented  grid.  Most  of  the  rotor  wake  lies  outside  of  such  a  local  grid  and  must 
be  accounted  for. 

The  first  treatment  of  hover  problems  with  CFD  methods  involved  the  use  of  a  small- 
disturbance  potential  equation.  The  representation  of  a  system  of  tip  vortices  passing  through 
the  grid  was  made  by  specifying  a  system  of  additional  constant-strength  branch-cut  sheets 
(see  Fig.  16)  using  the  same  logic  which  is  already  required  to  specify  the  rotor  trailing  sheet. 
The  edges  of  the  sheets  are  chosen  to  coincide  with  tip>- vortex  locations.  Note  that  in  Fig.  16 
the  additional  branch  cuts  are  shown  to  be  vertical.  Actually  the  sheet  orientation  is  irrelevant 
because  it  is  only  the  edge  of  a  constant-strength  cut  that  represents  the  vortex.  The  use  of 
branch  cuts  to  specify  shed  vortices  is  especially  easy  for  a  Cartesian  or  “H-type”  grid  (which 
all  small-disturbance  methods  use)  because  the  tip  vortices  are  nearly  parallel  to  the  grid 
lines.  This  approach  cannot  easily  be  used  to  specify  an  inboard  vorticity  sheet,  however, 
because  it  would  not  usually  coincide  with  a  coordinate  plane.  Fortunately  these  sheets 
are  weak  compared  to  the  tip  vortices  and  can  be  excluded  from  the  near-blade  problem. 
But  these  sheets  and  all  the  wake  elements  not  contained  in  the  grid  induce  much  of  the 
inflow  and  need  to  be  treated.  In  Ref.  [38]  a  small-perturbation  near-fleld  computation  was 
coupled  to  the  boundary-integral  code  HOVER  [18].  The  coupling  involved  a  modification 
of  the  HOVER  inflow  computation  (the  Biot-Savart  integral)  wherein  all  the  wake  elements 
contained  in  the  grid  were  excised  to  produce  a  “partial  inflow”  which  was  then  applied  to 
the  blade  surface  boundary  condition  as  a  “partial  angle  of  attack.”  This  resulting  blade-load 
distribution  provides  the  wake-circulation  distribution  required  by  HOVER.  HOVER  also  used 
an  experimentally  measured  prescribed  wake  geometry  as  an  input. 

Clearly  the  main  issue  in  hybrid  methods  is  the  manner  of  communicating  the  wake  data 
to  the  CFD  grid.  For  convenience,  in  future  discussions  we  will  refer  to  the  present  branch-cut 
approach  as  “gamma  coupling.”  The  use  of  a  surface  partial  inflow  is  referred  to  as  “alpha 
coupling.”  The  above  approach  used  both  methods. 

Although  gamma  coupling  is  simple  when  a  vortex  can  be  aligned  with  a  grid  line,  it  is 
not  easily  applied  to  arbitrary  grid/vortex  orientations.  One  possible  approach  would  be  to 
use  adaptive  grids;  but  there  are  simpler  ways  to  treat  the  problem. 

One  alternate  approach  involves  modifying  the  flow  equations.  We  can  represent  the  flow 
velocity  as 

V  =  Voo+Q  +  V0  (51) 

where  Voo  is  the  undisturbed  free-stream  velocity  that  results  from  translation  and  rotation, 
and  Q  is  an  induced  flow  field.  In  this  case  Q  is  the  velocity  induced  by  the  additional  wake 
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elements  in  the  grid.  Substituting  Eq.  (51)  into  the  continuity  equation  (Eq.  (1))  using  an 
inertial  frame  for  simplicity,  we  have 

Pt  +  V  •  (pW)  =  -V  •  pQ  (52) 

It  can  be  seen  that  the  only  effect  on  the  mass  equation  is  to  add  a  known  forcing  function 
to  the  right-hand  side.  Note  also  that  Q  can  be  a  rotational  field  with  no  loss  of  validity  of 
Eq.  (52).  However,  the  Bernoulli  equation  would  no  longer  be  valid  in  these  rotational  regions. 
But  these  regions  tend  to  be  thin  and  their  primary  importance  is  their  circulation.  Pressure  or 
density  variations  in  these  regions  have  no  global  effect.  This  velocity  decomposition  approach 
was  first  used  by  Steinhoff  [39,  40].  We  will  refer  to  the  use  of  an  induced  velocity  field  to 
specify  wake  systems  as  “Q  coupling.”  Two  different  Q-coupled  hybrid  hover  analyses  will 
now  be  described. 

Q  coupling  has  been  used  to  hybridize  the  FPR  and  HOVER  codes.  Because  FPR  uses 
a  spanwise  stacked  “O  grid,”  the  branch-cut  mode  of  vortex  representation  is  not  feasible. 
However,  the  velocity  decomposition  approach  (Q  coupling)  has  been  successfully  implemented 
[31].  In  this  analysis,  the  effect  of  the  far  wake  was  analyzed  using  a  partial  angle-of-attack 
approach  (alpha  coupling).  The  close  tip- vortex  elements  (those  passing  through  the  grid) 
were  approximated  by  infinite  straight  vortices,  which  then  defined  the  Q  field.  (Straight 
vortices  are  a  reasonable  approximation  and  obviate  a  Biot-Savart  computation).  Of  course 
Q  coupling  also  implies  its  own  surface  inflow  modification  because  changing  the  velocity  field 
entails  a  corresponding  change  in  the  no-flow-through  condition. 

Q  coupling  has  also  been  employed  by  Egolf  and  Sparks  [41]  using  the  nonconservative 
code  ROT22.  This  hover  implementation  differs  from  the  previous  analyses  in  that  the  effect 
of  the  total  wake  is  specified  on  the  entire  outer  computational  boundary,  which  eliminates 
the  need  for  a  partial  inflow  on  the  rotor  surface.  This  has  the  advantage  that  it  obviates  any 
need  to  modify  the  inflow  prediction  program.  Another  interesting  feature  of  this  work  is  that 
it  employed  an  adaptive  grid,  which  permitted  the  shed  wake  (from  the  trailing  edge  to  the 
rear  grid  boundary)  to  convect  freely.  This  work  also  uses  a  local  line  vortex  representation 
to  define  a  Q  field. 

The  above  discussion  presented  a  variety  of  finite-difference  rotor  computations  including 
an  extensive  self-induced  wake  system.  The  described  analyses  are  really  “pre-engineering” 
pilot  methods,  intended  to  demonstrate  the  feasibility  of  combining  difference  methods  with 
our  present  integral  methods;  and  they  require  validation  data  which  can  provide  both  wake 
geometry  (for  inflow  prediction)  and  high-speed  blade  surface  pressures. 

Toward  this  end,  an  extensively  pressure  instrumented  model  rotor  was  hover  tested 
[42]  by  the  (then)  U.S.  Army  Aeromechanics  Laboratory.  Simultaneous  measurements  of 
the  wake  geometry  (depicted  in  Fig.  17)  provides  the  tipn vortex-location  information  that  is 
needed  for  computational  vortex-modeling  studies.  This  data  has  been  studied  using  the  three 
previously  mentioned  hybrid  analyses,  i.e.,  the  small-perturbation  approach  which  combined 
gamma  and  alpha  coupling  [39],  and  conservative  (FPR  [31])  and  nonconservative  (ROT22 
[41])  full- potential  schemes  using  Q  coupling.  These  methods  all  assume  that  the  vortex 
strength  is  equal  to  the  maximum  blade  bound  circulation.  Figure  17  shows  a  comparison  of 
all  three  approaches,  with  the  data.  Overall,  there  is  a  remarkable  agreement  between  the  three 
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approaches.  In  all  three  methods  the  predicted  shock  location  is  slightly  aft  of  the  measured 
location.  This  discrepancy  is  probably  a  result  of  the  isentropic  flow  approximation  [38].  The 
agreement  between  the  conservative  and  nonconservative  predicted  shock  location  is  somewhat 
surprising  and  can  result  from  any  number  of  computational  details.  The  nonconservative 
approach  almost  always  produces  a  weaker  shock  than  its  conservative  counterpart  does.  The 
important  point  is  that  the  basic  soundness  of  the  blade/wake  aerodynamic  matching  process 
is  shown.  Accordingly,  it  appears  that  we  should  be  able  to  effectively  combine  our  flnite- 
difference  blade  computations  with  existing  boundary-integral  codes.  Such  analyses  could 
easily  include  fuselage-induced  upwash  effects.  Of  course,  in  order  for  the  analysis  to  be 
meaningful,  the  difference  portion  will  have  to  be  able  to  predict  drag. 

The  high  speeds  involved  in  the  previous  hover  example  are  not  operationally  realistic, 
and  transonic  flow  is  not  a  major  hover  consideration,  although  it  does  occur.  In  forward 
flight,  however,  such  tip  speeds  are  common  and  the  use  of  transonic  flow  models  becomes 
important.  In  principle,  the  use  of  a  near-fleld  wake  representation  (one  contained  in  the  grid) 
combined  with  a  far-fleld  model  should  permit  complete  forward-flight  analyses.  However,  the 
practical  difficulties  of  such  a  computation  are  much  greater  than  those  for  hover.  The  fact 
that  the  wake/blade  orientation  is  time- varying  appears  to  preclude  the  use  of  a  branch-cut 
wake  representation.  Furthermore,  the  use  of  velocity  decomposition  (Q  coupling)  requires  a 
Biot-Savart  integral  for  every  grid  point  (or  for  enough  grid  points  that  interpolation  can  be 
used)  and  for  every  time  step — a  rather  expensive  proposition  with  present  techniques.  Also 
the  straight- vortex  model  used  for  hover  is  probably  not  a  good  model  for  curved  vortices  in 
forward  flight.  For  these  reasons  a  full  forward-flight,  velocity-coupled  computation  has  not 
yet  been  performed,  although  it  is  obviously  possible. 

There  are  also  physical  reasons  to  defer  the  use  of  a  specified  complete  wake  model  in 
the  current  CFD  codes.  First,  the  wake  structure  of  even  high-advance-ratio  rotors  is  not  well 
known.  At  high  advance  ratios,  the  wake  can  be  weaker  (especially  on  the  advancing  side) 
and  farther  removed  from  the  blade.  Furthermore,  the  wake-induced  inflow  can  be  a  small 
percentage  of  the  total  inflow.  It  thus  becomes  reasonable  to  use  a  simpler  wake  model  for 
many  flight  conditions. 

At  the  present  time,  therefore,  forward-flight  hybrid  computations  have  been  performed 
only  on  the  basis  of  partial-angle-of- attack  or  alpha  coupling  (see  Fig.  18).  That  is,  the  wake 
elements  that  pass  through  the  computational  grid  are  not  modeled  by  means  of  their  spacial 
velocity  fields.  Instead,  their  inflow  is  accounted  for  only  in  the  body  boundary  condition  as 
a  partial  angle  of  attack.  This  angle  of  attack  is  still  “partial”  because  any  finite-difference 
computation  includes  a  circulation  sheet  emanating  from  the  trailing  edge.  The  portion  of  the 
wake  that  corresponds  to  the  grid  branch-cut  region  must  be  excised  from  the  Biot-Savart 
integration  in  the  outer  boundary-integral  wake  model.  Failure  to  perform  this  modification 
would  result  in  a  double  accounting  for  the  shed  wake  circulation. 

In  computing  these  advancing  blade/ wake  flows  it  is  convenient  to  include  all  geometric 
(twist),  blade  motion  (flapping  and  deformation),  and  inflow  effects  in  this  partial  angle  of 
attack.  A  fixed  untwisted  grid  is  commonly  used  and  the  partial  angle  of  attack  is  specified  by 
applying  a  flow-through  (transpiration)  boundary  condition  on  the  surface.  The  alternative 
(for  a  code  with  a  body-conforming  grid)  would  be  to  generate  a  new  grid  (at  each  time  step) 
that  included  the  partial  angle  as  an  effective  twist.  On  this  grid,  a  no-flow-through  boundary 
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condition  would  then  be  used.  This  approach  is  much  more  difficult,  however,  and  does  not 
significantly  increase  accuracy. 

Up  to  this  point  we  have  referred  to  hybridizing  only  in  relation  to  different  aerodynamic 
models.  However,  we  can  only  confine  our  discussion  to  aerodynamics  when  we  know  the 
boundary  conditions,  and  for  an  advancing  rotor  we  do  not  know  the  blade  motion  beforehand. 
Instead,  we  typically  know  the  primary  (lift  and  propulsive)  forces  and  must  then  perform  an 
iterative  aerodynamic/dynamic/elastic  computation  to  find  the  control  inputs  (e.g.,  cyclic  and 
collective)  that  produce  these  forces.  This  is  referred  to  as  the  rotor  trim  problem.  The  end 
result  of  the  trim  process  is  a  knowledge  of  the  blade  motion.  A  typical  trim  process  will  use 
2-D  airfoil  table  look-ups  to  provide  the  local  blade  aerodynamics  around  the  entire  gizimuth. 
Such  a  process  can  require  the  computation  of  about  10  rotor  revolutions  to  converge  to  a 
solution.  This  is  a  fast  process  with  look-up  tables,  but  it  is  clearly  impractical  to  directly 
replace  the  look-up  process  by  CFD  solutions  at  each  time  step.  The  convergence  of  present 
trim  procedures  is  such  that  CFD  solutions  cannot  be  placed  in  the  trim  loop. 

A  solution  to  this  problem  that  still  uses  present  tabular  methods  is  to  solve  the  local 
rotor  problem  outside  of  the  trim  loop.  This  difference  solution  is  then  used  as  a  brse  about 
which  to  find  corrections  resulting  from  the  trim  process.  Therefore  the  lift  is  computed  as 
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where  a  and  Qoid  are  the  angles  of  attack  from  the  current  and  previous  trim  loops,  respectively. 
The  solution  has  converged  when  a  — ♦  aoid  and  the  lift  correction  vanishes.  At  this  point  the 
finite-difference  computed  lift  is  fully  consistent  with  the  rotor  inflow  and  motion.  This  scheme 
appears  at  first  to  be  a  slight  retreat  from  the  goal  of  obviating  table  look-ups.  However,  the 
tables  affect  only  the  convergence  rate,  not  the  final  answer. 

This  forward-flight  matching  process  was  first  performed  [43]  using  a  small-disturbance 
code  (FDR)  and  the  comprehensive  code  CAMRAD.  Currently,  the  same  matching  is  regularly 
performed  using  the  conservative  full-potential  code  FPR.  The  matching  of  the  CAMRAD 
and  FPR  codes  is  summarized  in  Fig.  19.  The  process  is  started  by  obtaining  a  trimmed 
nonuniform  inflow  solution  (with  a  full  vortex-lattice-modeled  wake)  with  the  lift  obtained 
totally  from  airfoil  tables.  This  is  the  normal  operation  of  CAMRAD  except  that  partial 
angles  of  attack  are  also  computed.  These  partial  angles  provide  the  necessary  boundary 
conditions  to  obtain  an  FPR  solution  for  the  lift  distribution.  The  program  then  alternately 
computes  new  modified  trim  solutions  and  FPR  solutions,  until  the  lift  correction  vanishes. 
This  is  an  efficient  scheme  since  it  does  not  do  the  most  time-consuming  tasks  (influence 
coefficients  and  finite-difference  computations)  in  the  innermost  trim  loop.  The  convergence 
of  this  process  is  very  fast.  Stiff  rectangular  rotors  frequently  give  good  results  in  one  iteration. 
Soft  rotors  with  varying  planforms  have  required  about  three  iterations. 

An  example  of  such  a  computation  is  illustrated  in  Fig.  20,  which  shows  a  comparison  with 
one  of  the  ONERA  three-bladed  test  cases.  For  this  case,  the  tip  rotational  Mach  number 
is  0.628  and  the  advance  ratio  (ratio  of  forward  speed  to  rotational  spi  J)  is  0.387.  The 
airfoils  used  were  S-130XX  (a  variant  of  the  NACA  five-digit  family).  This  computa*’on  was 
accomplished  using  the  coupled  FPR  and  CAMRAD  codes  [44].  Overall,  the  agreement  of 
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the  data  and  the  computation  is  good.  However,  the  computation  somewhat  overpredicts 
the  upper  pressures  early  in  the  first  quadrant  (0°  <  V’  <  90°)  and  the  lower  surfaces  in  the 
second  quadrant.  The  shock  location  and  strength  seem  to  be  well-predicted  throughout  the 
computation. 

The  previous  cases  involved  stiff,  articulated  model  rotors,  in  which  the  blade  motion  (ex¬ 
cept  for  rotation)  was  almost  pure  flapping  with  very  little  elastic  torsion  or  bending.  Full-scale 
rotors  tend  to  be  much  softer  than  the  above-mentioned  model  blades,  and  elastic  deforma¬ 
tions  can  be  important.  The  first  published  computational/experimental  data  comparison  for 
an  actual  flight  vehicle  involved  an  Aerospatiale  SA349.  A  typical  data/computation  compar¬ 
ison  [45],  performed  using  CAMRAD/FPR,  is  shown  in  Fig.  21.  Although  this  comparison 
of  surface  pressures  is  promising,  the  differences  between  the  computed  and  measured  lift 
actually  constitute  a  significant  load  error.  The  resolutions  of  these  problems  will  probably 
require  improvements  in  our  structural,  wake,  and  coupling  models.  Nevertheless,  the  present 
coupling  has  effectively  integrated  3-D,  unsteady  transonic  flow  analysis  into  the  comprehen¬ 
sive  rotor  modeling  process.  These  codes  constitute  a  valuable  analysis-and-advanced-design 
tool  in  the  hands  of  a  knowledgeable  user. 

One  of  the  areas  that  requires  further  review  is  the  basic  notion  of  the  partial  angle-of- 
attack  coupling.  The  idea  of  a  chordwise  constant  inflow  is  only  valid  if  all  the  wake  elements 
causing  that  inflow  are  well-removed  from  the  blade  (i.e.,  by  more  than  a  chord).  Recent 
studies  of  blade/vortex  interactions  (nearly  direct  impingement)  have  shown  good  correlations 
with  BVI  leading-edge  pressure  data  merely  by  using  a  chordwise- varying  inflow  [46].  Other 
computations  show  significant  differences  between  the  use  of  a  chordwise-varying  inflow  and 
a  full  Q-field  representation  of  the  vortex  [47].  (See  the  related  BVI  discussion  in  Section 
V.)  Perhaps  the  most  significant  element  of  these  studies  concerns  the  effects  of  unsteadiness. 
The  present  angle-of-attack  (or  chordwise-constant  inflow)  coupling  lumps  all  unsteady  effects 
(except  for  Mach-number  variation)  into  a  single  surface  inflow.  No  diflferentiation  is  made 
between  inflow  and  flap/pitch  variation.  It  is  well  known  that  the  response  to  these  two  types 
of  excitation  is  not  the  same.  Recent  computational  correlation  studies  with  in-flight  data 
[48]  indicate  stronger  discrepancies  in  the  unsteady-flow  model  than  previously  expected.  Of 
course,  such  discussions  would  be  moot  if  we  had  a  unified  representation  of  the  rotor /wake 
system. 

UNIFIED  FLOW  METHODS 

The  rotor  and  its  wake  constitute  a  tightly  knit  system,  and  it  is  natural  to  solve  it  as 
a  single  problem.  This  has  been  done  for  years  using  boundary-integral  methods.  The  need 
to  include  compressibility  has  driven  efforts  to  do  the  same  with  differential  methods,  using 
a  single  grid  that  encompasses  both  the  blade  and  the  wake.  This  work  has  excluded  con¬ 
ventional  potential  methods  because  of  an  assumed  inability  to  convect  the  wake  freely.  Thus 
potential  methods  had  been  relegated  to  use  as  local-blade-flow  solvers  in  hybrid  analyses, 
and  there  have  been  various  attempts  to  use  conventional  Euler  solvers  for  unified  CFD  rotor 
modeling. 

Euler  solvers  have  not  been  discussed  to  this  point.  It  suffices  for  present  purposes 
to  mention  that  many  such  solvers  now  exist.  Most  use  the  finite- volume  approach  and  an 
explicit,  centered-difference  solution  scheme,  which  employs  user-specified  dissipation  terms  to 
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ensure  stability  [49].  Although  these  methods  require  “tuning”  of  dissipation  terms  (required 
for  stability),  they  are  the  most  efficient  Euler  solvers.  Another  class  of  Euler  solver  uses 
implicit,  upwind-differenced  schemes  that  require  no  tuning  and  appear  to  be  less  dissipative. 
Several  of  these  Euler  methods  have  also  been  used  in  hybrid  computations,  in  the  manner 
described  in  the  previous  section  [50]. 

The  first  unified  blade/wake  Euler  computation  of  a  hovering  rotor  was  performed  by 
Kroll  [51]  using  a  centered  explicit  solver.  This  computation  modeled  the  hover  experiment  of 
Ref.  [42]  using  an  “0-0”  grid  fixed  to  a  blade  and  large  enough  to  encompass  a  considerable 
portion  of  the  wake.  A  problem  found  with  this  computation  was  that  the  wake  was  dissipated 
to  such  an  extent  that  the  induced  inflow  was  underpredicted.  That  dissipation  was  a  problem 
was  shown  by  performing  a  grid  sensitivity  study.  With  cr  arse  grids  the  lift  distribution  was 
badly  overpredicted,  especially  inboard,  indicating  an  underprediction  of  inflow.  Refining  the 
grid  reduced  the  lift,  but  not  sufficiently  to  compare  well  with  the  experimental  data  (Fig.  22). 
A  comparison  of  the  computed  pressure  distributions  with  the  data  was  very  favorable  at  the 
tip,  however.  Similar  efforts  by  Kramer  et  al.  [52],  and  Chen  and  McCroskey  [53]  used  an 
upwinded  Euler  scheme  to  make  comparisons  with  the  same  hover  data,  with  similar  results 
(very  good  agreement  with  the  outboard  data,  but  a  tendency  to  overpredict  the  lift  inboard). 
Chen’s  results  showed  only  a  slight  grid-refinement  sensitivity.  This  was  probably  a  result  of 
the  far  grid  (away  from  the  blade)  being  so  stretched  (to  avoid  excessive  computing  time)  that 
the  vortex  was  still  unavoidably  dissipated.  This  vortex  dissipation  can  be  seen  in  Fig.  23, 
which  shows  computed  vorticity  contours  at  various  distances  from  the  trailing  edge.  With 
present  algorithms,  the  cost  of  such  Euler  computations  is  very  high,  and  it  is  prohibitive 
with  a  grid  which  is  dense  enough  to  eliminate  the  dissipation.  The  use  of  unified  rotor/ wake 
computations  thus  seems  remote,  but  new  developments  using  potential  methods  as  part  of 
the  solver  could  change  this  situation. 

Steinhoff  and  Ramachandran  [54,  55,  56]  have  recently  re-examined  the  potential  ap¬ 
proach  and  arrived  at  a  new  method  of  specifying  the  shed  circulation  such  that  it  becomes 
free  to  convect  through  the  grid.  The  heart  of  their  approach  is  the  use  of  velocity  decompo¬ 
sition  (Eq.  (51))  which  was  cited  in  connection  with  the  earlier  discussion  of  hybrid  methods. 
The  central  idea  of  this  work  is  that  when  a  Q  field  is  used  to  represent  a  wake,  any  form  of 
velocity  field  Q  can  be  used  as  long  as  V  x  Q  =  uf,  where  uf  is  the  vorticity  distribution  for 
the  wake  sheet.  This  new  formulation  thus  entails  spreading  the  circulation  sheet  to  give  it 
an  at  distribution,  and  then  finding  an  appropriate  Q,  which  defines  <j.  This  Q  then  defines  a 
forcing  function  for  the  mass  equation  (Eq.  (48)).  This  process  has  been  previously  described, 
but  in  this  work  a  different  form  is  chosen  for  Q.  Using  the  Biot-Savart  law  would  allow  us  to 
find  Q  as  the  velocity  induced  by  u.  This  is  a  field  that  fills  all  space  and  is  therefore  too  costly 
to  compute,  since  each  field  point  would  require  an  integration  over  the  entire  sheet.  However, 
Q  can  also  be  defined  as  a  velocity  normal  to  the  sheet  that  defines  u)  (see  Fig.  24).  Such 
a  field  is  zero  everywhere,  except  where  ui  is  nonzero  and  therefore  is  computed  at  relatively 
few  points.  This  representation  seems  unphysical  at  first,  because  such  normal  velocities  do 
not  exist  in  the  wake.  However,  on  finding  (f)  (from  the  continuity  equation)  and  adding  all 
velocity  components  to  get  the  total  V,  this  normal  velocity  is  cancelled  by  V<f),  leaving  the 
expected  induced  flow  throughout  the  entire  field.  That  this  must  happen  becomes  obvious 
when  we  recall  that  the  Biot-Savart  law  is  itself  derived  from  the  continuity  equation  (i.e.. 
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Laplace’s  equation).  Prom  this  point  of  view  we  could  think  of  this  potential  solver  as  a  non¬ 
linear,  Biot-Savart  induced-flow  solver.  Moreover,  this  is  the  same  solver  that  is  solving  for 
the  local  flow  on  the  rotor.  Subsequently,  the  symbol  will  be  used  to  denote  this  normal 
form  of  Q. 

To  And  the  required  strength  of  we  use  Gauss’s  theorem  to  obtain  a  relation  from  the 
integral  of  along  a  normal  through  each  point  on  the  sheet,  thus 


r  = 


(54) 


The  circulation,  P,  is  known  at  the  upstream  edge  of  the  sheet  (blade  trailing  edge)  from 
the  lift  distribution,  which  is  computed  as  part  of  the  entire  calculation.  Since  T  is  constant 
along  mean  streamlines  within  the  sheet,  it  can  easily  be  computed  on  the  entire  sheet.  This 
r  distribution  provides  a  scaling  factor  which  gives  the  magnitude  of  as  soon  as  the  width 
and  functional  form  of  Q'^  are  determined.  This  width  and  functional  form  can  be  found  by 
means  of  a  viscous  solution  or  by  simply  choosing  computationally  convenient  forms.  For  the 
hover  wake  problem,  the  latter  approach  suffices.  A  particularly  useful  form  for  involves  a 
Clebsch-type  [57]  representation; 

Q"  =  rVA  (55) 

where  is  a  3-D  field  which  smoothly  goes  to  the  appropriate  P  on  the  sheet  as  r  approaches 
the  sheet  surface.  In  this  representation,  P  is  a  local-sheet-strength  function  (which  is  actually 
smeared  spacially  to  facilitate  the  treatment  of  highly  curved  sheet  surfaces)  and  VA  is  a  spacial 
distribution  function.  The  term  A  is  any  convenient  function  whose  gradient  is  nonzero  within 
some  specified  small  distance  from  the  sheet.  This  specified  smearing  distance  is  chosen  to  be 
on  the  order  of  the  local  grid  size. 

The  major  computational  work  in  this  method  is  the  solution  of  the  continuity  equation 
for  the  potential.  To  solve  for  0,  a  finite-volume,  conservative,  full-potential  solver  is  used 
that  is  semi-implicit  (i.e.,  explicit  in  the  radial  direction  only)  and  employs  an  AF  scheme  in 
each  radial  plane.  A  blade-fixed  H  grid  is  used  to  solve  for  the  potential  (see  [55]  for  details). 
The  resulting  solution  involves  iteration  between  solving  the  mass  equation  and  convecting 
the  shed  circulation.  Since  this  is  a  Lagrangian  convection  process  there  is  no  possibility  of 
vorticity  dissipation.  During  the  computation,  a  four-step  procedure  is  repeatedly  used:  (1) 
the  vortex-sheet  position  is  integrated  as  a  set  of  marker  streamlines  to  follow  the  flow  using 
interpolated  values  of  V  from  the  fixed  grid;  (2)  Q'’  is  computed  at  grid  points  near  the  sheet; 
(3)  a  potential  0  is  computed  at  all  grid  points  by  solving  the  compressible  mass  conservation 
equation;  and  (4)  a  new  velocity  V  is  computed  at  each  grid  point  after  adding  Q'’  to  the 
potential  gradient  and  free-stream  components  of  the  velocity.  At  convergence  the  vortex 
sheet  follows  the  flow.  This  procedure  has  been  implemented  in  a  code  called  HELIX-I. 

This  new  solution  procedure  has  been  applied  to  a  number  of  hovering  rotor  configura¬ 
tions.  Comparisons  have  been  made  with  experimental  measurements  of  tip>-vortex  geometry, 
thrust,  and  power.  In  order  to  compute  the  latter  it  was  necessary  to  use  a  simple  integral 
boundary-layer  code  for  skin  friction  and  an  energy-flux  integral  for  induced  power.  Figure 
25  shows  the  computed  wake  geometry  for  a  4-bladed,  linearly  twisted  model  rotor  [13].  The 
wake  is  seen  both  by  the  wake  marker  loci  (at  a  point  25°  behind  a  blade)  and  the  contours  of 
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vorticity  that  these  markers  carry  through  the  field.  The  tip  vortices  can  be  clearly  seen.  The 
axial  and  radial  convection  rates  of  these  vortices  are  compared  with  experiment  in  Figs.  26 
and  27  and  the  agreement  is  good.  The  experimental-to-computational  comparison  of  the 
lift/power  polars  is  shown  in  Fig.  28.  The  computations  somewhat  overpredict  the  power  at 
higher  lift,  for  reasons  that  are  not  yet  clear.  These  are  the  first  performance  polars  ever 
predicted  by  differential  CFD  methods.  Figures  29  and  30  show  a  similar  set  of  comparisons 
for  a  model  Boeing  360  rotor,  and  the  comparisons  are  favorable  [58]. 

This  is  a  new  approach  to  flow  computation,  which  combines  the  traditional  CFD  ability 
to  treat  supercritical  flows  with  the  Lagrangian  convection  typical  of  the  integral  codes.  Fur¬ 
thermore,  the  approach  greatly  expands  the  capabilities  of  potential  methods.  The  approach 
is  extendable  to  forward  flight  and  should  be  able  to  handle  any  of  the  previously  mentioned 
advancing-rotor  problems,  but  development  work  is  required.  It  is  significant  that  this  class  of 
methods  is  the  first  in  which  differential  CFD  is  used  to  produce  hover  results  of  engineering 
interest. 

V.  ADVANCED  ROTOR  FLOW  TOPICS 

The  previous  section  dealt  with  integrating  the  local  blade  analysis  with  the  wake.  With 
the  success  of  this  integration  in  view,  we  now  review  some  previously  excluded  flow  topics 
that  need  to  be  included  in  the  total  rotorcraft  model  in  order  to  convert  CFD  into  a  truely 
complete  analysis. 

DYNAMIC  STALL  AND  OTHER  SEPARATED  FLOWS 

Retreating  blade  dynamic  stall  is  one  of  the  most  dangerous  rotor  flow  conditions.  It  is 
also  one  of  the  oldest  aerodynamic  problems  involving  major  computational  effort. 

The  oldest  line  of  computational  stall  modeling  used  boundary-element  analyses  to  model 
the  unsteady,  separated  flow  on  dynamically  stalling  airfoils  [59,  60].  Such  codes  represent 
the  profile  given  by  a  vortex  lattice  or  panel  model  and  permit  a  circulation  sheet  to  peel 
from  the  blade  surface,  on  cue  from  a  boundary-layer  analysis.  The  sheet  organizes  itself  into 
something  similar  to  a  stall  vortex,  convects  downstream,  and  produces  loads  that  resemble 
measured  values.  The  main  problem  is  the  accurate  determination  of  the  stall  commencement 
point  and  time. 

One  of  the  first  stall  analyses  in  the  modern  CFD  sense  is  that  of  Wu  et  al.  [61,  62].  This 
approach  is  a  combined  integral-differential  scheme  in  which  the  vorticity  diffusion  equation 
is  solved  by  difference  methods  and  is  simultaneously  coupled  to  a  Biot-Savart  integral  for 
the  induced  velocity  field.  This  approach  can  be  especially  efficient  because  the  formulation 
only  requires  computation  where  vorticity  is  nonzero.  Figure  31,  which  illustrates  the  result 
from  such  a  computation,  shows  the  streamlines  and  vorticity  contours  for  a  NACA0012 
airfoil  undergoing  sinusoidal  oscillations;  a  =  15°  -t-  10sin(u;t),  where  reduced  frequency  k  = 
ujc/2V  =  0.15.  These  computations  clearly  show  the  separation  and  downstream  convection 
of  the  classical  leading-edge  stall  vortex.  Flow  reattachment  can  be  seen  also.  These  events 
are  reflected  in  the  lift  and  moment  plots  of  fig.  32.  Unfortunately,  the  Biot-Savart  integral 
limits  this  approach  to  incompressible  flows.  This  is  a  serious  limitation,  even  for  the  low 
Mach  numbers  which  characterize  retreating  blade  stall  (Mach  0.3-0. 5). 
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One  of  the  first  purely  differential  CFD  stall  analyses  was  that  of  Mehta  [63],  who  per¬ 
formed  laminar,  incompressible  computations  that  produced  good  qualitative  comparisons 
with  stall  data  (see  Fig.  33). 

Recently,  a  more  general  compressible  Reynolds-averaged,  Navier-Stokes  CFD  method 
has  been  developed  by  Sankar  et  al.  [64].  This  involves  a  centered,  finite-difference,  implicit 
scheme  that  treats  the  viscous  terms  explicitly.  Like  all  centered  schemes,  the  method  uses 
specified  artificial  dissipation  terms  to  assure  stability  (see  [65]).  Figure  34  from  [66]  shows  a 
comparison  of  computed  and  measured  lift  and  pitching  moment  for  a  sinusoidally  oscillating, 

2- D  airfoil  (Sikorski  SSC-A09).  For  this  case.  Moo  =  0-2,  Re  =  2.0 x  10®,  and  k  =  ujcf2V  =  0.1. 
The  angle  of  attack  varies  between  0°  and  20°.  Generally  good  agreement  is  obtained  for  both 
lift  and  moment,  but  there  are  important  differences.  The  computed  stall  conunences  sooner 
than  the  measured  stall,  which  begins  at  the  top  of  the  stroke.  The  earlier  computed  stall  can 
be  seen  in  both  the  pitching-moment  drop  and  in  the  oscillations  in  the  lift.  The  peak  loads 
are  well  predicted,  especially  the  maximum  nose-down  moment.  However,  this  is  a  deep-stall 
case  (one  dominated  by  a  distinct  vortex  generated  at  the  leading  edge),  and  many  methods 
that  generate  a  leading-edge  vortex  (including  a  vortex-lattice  method)  will  reproduce  similar 
peak  moments.  The  problem  is  to  accurately  predict  the  beginning  of  stall.  This  becomes 
especially  important  for  light  stalls,  in  which  the  well-organized  vortex  does  not  occur.  Such 
stalls  can  generate  the  greatest  negative  pitch  damping  and,  hence,  stall  flutter.  All  published 
computational-experimental  comparisons  (that  this  writer  is  aware  of)  treat  only  deep  stall.  A 
systematic  comparison  of  loads  and  pitch  damping  for  light  stalls  would  be  useful.  This  could 
be  difficult  because  of  the  present  state  of  turbulence  modeling.  The  above  computation  used  a 
Baldwin-Lomax  (BL)  algebraic  turbulence  model  that  is  not  intended  for  strongly  separated 
flows.  In  [64]  a  K-e  model  was  compared  to  a  BL  model  for  deep  stall  computations,  and 
surprisingly  little  improvement  was  found.  It  seems,  then,  that  the  basic  impediments  to  the 
solution  of  this  problem  are  physical  more  than  numerical.  Even  our  ability  to  predict  steady 
stall  is  not  well  demonstrated.  The  examples  shown  above  are  a  major  advance,  but  this 
advance  is  mainly  a  result  of  our  improving  computational  capabilities. 

Our  computational  ability  to  model  unsteady  stall  does  enable  us  to  predict  some  major 
trends  that  are  not  always  or  entirely  governed  by  the  details  of  turbulence.  For  example. 
Figs.  35  and  36  show  the  computed  and  experimental  lift  history  for  an  airfoil  (SSC-A09) 
undergoing  a  constant  rate  pitch- up  [66].  It  can  be  seen  in  these  cases  that  the  effects  of  Mach 
number  and  pitch  rate  on  stall  initiation  are  well  predicted,  although  the  subsequent  events 
are  not.  It  might  be  useful  to  seek  other  effects  whose  trends  can  be  predicted  by  present 
computational  methods.  For  example,  it  is  well  known  that  most  stall  processes  involve  large 

3- D  effects.  Furthermore,  major  claims  have  recently  been  made  concerning  the  usefulness 
of  nonrectangular  planforms  for  the  alleviation  of  stall  effects  on  very-high-speed  helicopters. 
An  ability  to  demonstrate  this  usefulness,  both  computationally  and  experimentally,  would  be 
of  major  scientific  and  engineering  importance.  To  date,  only  steady  computations  [67]  have 
been  maxle  for  such  configurations.  Figure  37  shows  a  numerical  visualization  of  the  computed 
flow  on  a  nonrotating,  nonrectanglar  blade  using  simulated  particle  trajectories.  This  TNS 
computation  shows  tip-vortex  formation,  inboard  separation  regions,  and  general  qualitative 
agreement  with  observed  flows  (oil-flow  visualizations). 
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Steady  3-D  TNS  computations  such  as  the  preceding  can  also  be  useful  in  their  own 
right,  especially  for  illuminating  the  details  of  tip-vortex  formation.  The  tip  vortex  interacts 
with  following  rotor  blades  and  is  therefore  an  important  element  in  the  prediction  of  higher- 
harmonic  loads  and  noise.  Furthermore,  it  has  been  shown  that  the  vortex  formation  process 
can  have  an  important  influence  on  the  tip  drag  [68,  69].  This  process  has  been  studied  nu¬ 
merically  by  Srinivasan  et  a!.  [70],  again  using  the  TNS  equations  with  an  algebraic  turbulence 
model.  The  effects  of  different  types  of  tip  cutoffs  on  the  vortex  formation  were  studied  and 
compared  to  experiment.  Figure  38  shows  a  numerical  flow- visualization  comparison  of  the  tip 
vortex  formation  for  a  rounded  and  a  flat  wingtip.  It  can  be  seen  that  the  sharp  edges  of  the 
squared  tip  induce  an  early  lift-off  of  the  tip  vortex.  Numerically  generated  patterns  of  surface 
particle  flows  have  been  compared  to  actual  oil-flow  visualizations,  as  shown  in  Fig.  39  for  a 
rounded  tip.  Figure  40  shows  a  comparison  of  computed  and  experimental  chordwise  pressure 
distributions  near  a  wingtip.  The  test  data  are  for  a  flat  tip  only,  whereas  the  computational 
results  are  for  flat,  round,  and  beveled  tips.  The  flat-tip  computations  have  all  the  qualitative 
features  of  the  data  and  become  quantitatively  correct  away  from  the  tip.  In  spite  of  the  dif¬ 
ferences  seen  in  the  tip  pressures,  the  computational  values  for  lift,  drag,  and  pitching  moment 
are  very  close  to  the  measured  quantities.  Clearly,  computational  tools  do  exist  for  studying 
these  effects,  but  the  fine  details  of  these  flows  can  not  yet  be  accurately  predicted.  This  may 
be  due  to  the  well-known  turbulence  modeling  deficiencies.  The  same  problems  that  are  seen 
in  viscous  tip  computations  will  probably  also  occur  in  efforts  to  predict  fuselage  flows.  These 
flows  must  be  simulated  to  obtain  downloads  in  hover,  and  drag  in  forward  flight.  That  these 
problems  are  of  great  importance  is  illustrated  by  the  major  performance  improvements  that 
can  result  from  drag- reduction  programs.  CFD,  together  with  major  improvements  in  our 
physical  flow  models,  can  be  an  important  part  of  such  design  improvement  work. 

BLADE- VORTEX  INTERACTIONS 

Rotors  encounter  vortices  under  a  wide  variety  of  circumstances  and  the  resulting  in¬ 
teractions  are  a  fundamental  source  of  vibratory  loading  and  noise.  The  problem  that  most 
simply  embodies  such  interactions  is  that  of  a  vortex  convecting  past  a  2-D  airfoil.  The  idea 
in  solving  such  a  problem  is  that  the  techniques  developed  should  be  directly  transferable  to 
3-D  problems. 

Srinivasan  et  al.  [70]  have  studied  the  2-D  blade/vortex  interaction  (BVI)  using  both 
small-disturbance  and  TNS  flow  models.  Using  the  TNS  equations,  an  attempt  was  also  made 
to  directly  convect  the  vortex  as  part  of  the  total  flow  field.  This  attempt  suffered  from  the 
expected  numerical  diffusion.  However,  successful  nondiffusing  computations  were  then  made 
using  the  Steinhoff  velocity  decomposition  method.  Figure  41  shows  computed  grids  (these  are 
adaptive),  surface-pressure  distributions,  and  Mach  contours  for  such  a  computed  close  BVI. 
This  is  a  high-Mach-number  case  in  which  the  vortex  passes  through  the  supersonic  region 
and  momentarily  bifurcates  the  shock.  Favorable  comparisons  were  also  made  with  pressure 
data  from  a  rotor  BVI  test.  Probably  the  most  innovative  CFD  treatment  of  the  2-D  BVI 
problem  is  that  of  Rai  [71]  who  developed  a  high-order,  upwind  TNS  scheme  whose  dissipation 
is  low  enough  to  permit  a  full  velocity  field  convection  of  the  vortex.  This  is  also  the  first 
method  capable  of  predicting  head-on  BVIs.  Figure  42  shows  the  grid  that  was  used  for  these 
computations.  It  includes  a  high-density  upstream  grid  region  which  minimizes  the  vortex 
dissipation  prior  to  the  BVI.  Figure  43  shows  the  convection  of  the  vortex  before,  during,  and 
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after  the  BVI.  The  post-BVI  vortex  has  a  duplex  structure,  which  has  been  seen  in  rotor  BVT 
testing  with  smoke  visualization.  Figure  44  shows  the  corresponding  pressure  contours.  These 
are  especially  significant  because  they  reveal  the  initiation  of  an  acoustic  wave  generated  by 
the  BVI.  This  method  is  clearly  a  valuable  tool  for  the  study  of  basic  rotor  acoustic  signal 
generation,  and  has  been  used  in  the  work  of  Baeder  [72].  Figure  45  shows  the  generation 
and  early  development  of  a  BVI-generated  acoustic  wave  for  a  supercritical  flow.  Rai  has  also 
made  comparisons  with  the  rotor  BVI  surface-pressure  data  of  [73].  Figures  46  and  47  include 
a  schematic  of  this  experiment  and  a  comparison  of  the  pressure  data  and  computations, 
which  shows  good  agreement.  The  computed  results  shown  here  agree  very  closely  with  the 
results  of  Srinivasan  et  al.  [70].  The  cost  of  this  upwind  scheme  probably  excludes  it  from 
consideration  as  a  practical  computation  method,  but  it  constitutes  a  standard  against  which 
other  methods  should  be  compared. 

The  first  3-D  BVI  computations  were  performed  by  Caradonna  et  al.  [46],  who  used  the 
FPR  code  with  a  velocity-decomposition  model  of  the  vortex  flow  field.  Good  agreement 
was  obtained  with  the  BVI  pressure  data  of  [74],  for  both  parallel  and  oblique  blade- vortex 
interaction  angles.  Figure  48  shows  a  typical  comparison  of  computed  and  measured  surface 
pressures  for  a  transonic,  oblique  BVI.  A  surprising  feature  of  this  work  is  that  good  agree¬ 
ment  was  seen  in  surface-pressure  comparisons  even  for  head-on  parallel  BVIs.  Because  the 
present  velocity  decomposition  requires  a  specified  vortex  structure,  there  is  no  mechanism 
allowing  for  the  alteration  of  this  structure.  Furthermore,  the  Bernoulli  equation  is  not  valid 
within  the  vortex.  Even  so,  good  surface  pressure  agreement  was  obtained.  This  indicates  the 
possibility  that  the  detailed  momentum  equations  may  not  have  to  be  solved  strictly  within 
the  vortex  for  some  types  of  vortical  interactions — a  potentially  significant  possibility  when 
these  methods  are  combined  with  the  potential,  vorticity-convection  methods,  such  as  those 
described  in  Section  IV.  Another  insight  gained  from  this  work  involves  the  relative  effective¬ 
ness  of  Q  coupling  and  alpha  coupling.  Both  Q-field  and  surface-inflow  methods  of  vortex 
specification  were  tested.  Although  the  surface  inflow  comparisons  with  data  showed  good 
agreement,  the  Q-field  method  was  consistently  better.  The  Q-field  method  was  found  to  be 
superior  for  parallel  and  oblique  interactions,  and  for  subcritical  and  supercritical  flows.  The 
reason  for  this  is  probably  that  a  surface-inflow  method  involves  approximating  the  effect  of 
the  vortex  as  an  effective  body  motion.  It  is  well  known  from  linear  theory  (for  instance, 
comparing  the  Theodorsen  and  Sears  functions)  that  the  response  of  a  wing  to  a  plunge  and 
a  gust  are  initially  quite  different.  The  above  comparison  (between  Q  coupling  and  alpha 
coupling)  is  directly  analogous.  It  appears,  therefore,  that  we  will  have  to  find  a  convenient 
way  to  Q-couple  our  rotor  blade  computations. 

The  computational  and  experimental  study  of  blade/vortex  interactions  will  probably 
continue  for  some  time,  because  it  is  providing  insight  into  important  physical  problems  and 
is  a  rich  source  of  directions  for  computational  modeling. 

INGREDIENTS  FOR  PERFORMANCE  PREDICTION 

In  order  to  reap  the  benefits  of  constructing  a  rotor/wake/fuselage  CFD  analysis,  it  is 
necessary  to  study  the  accuracy  of  drag  prediction.  With  CFD,  our  first  means  of  drag 
prediction  is  to  directly  integrate  the  surface- pressure  distributions.  Drag  prediction  is  a  very 
rigorous  test  of  a  numerical  method  because  drag  is  almost  always  a  small  number,  involving 
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cancellations  of  the  forces  from  different  parts  of  the  blade.  Conversely,  lift  is  relatively  easy  to 
predict.  The  accurate  prediction  of  drag  requires  careful  attention  both  to  numerical  accuracy 
and  to  the  physics.  Numerical  accuracy  is  usually  a  matter  of  assuring  that  the  differencing 
of  geometric  and  flow  quantities  are  consistent  with  each  other  [75]  and  that  “good”  grids  are 
used.  Whereas  consistency  is  something  that  is  supposedly  built  into  a  code,  grid  construction 
is  usually  a  matter  for  the  user  and  tends  to  be  an  art.  There  are  only  a  few  grid-construction 
rules:  (1)  grids  should  be  as  orthogonal  as  possible,  (2)  rapid  changes  in  grid  size  should  be 
avoided,  and  (3)  the  ratio  of  cell  lengths  should  be  as  close  as  possible  to  one.  Although 
there  are  a  number  of  grid  generators,  few  (or  none)  are  right  for  all  situations,  and  user 
intervention  is  the  norm.  Such  intervention  is  usually  required  because  of  geometric  or  grid- 
size  limitations.  The  end  result  is  often  a  less-than-ideal  compromise,  which  can  result  in 
very  slow  convergence  or,  worse,  inaccuracy.  An  acid  test  of  a  solution  is  its  drag  accuracy. 
Therefore  it  is  often  necessary  to  perform  nonlifting,  subcritical  test  computations,  because 
these  are  the  only  cases  for  which  we  know  the  theoretical  drag.  Physical  accuracy  is  a  matter 
of  having  a  proper  viscous  model.  For  most  rotor  flows,  a  noninteracting  boundary-layer 
model  is  probably  adequate.  For  supercritical  flows,  the  shock  entropy  rise  is  a  greater  and 
more  performance-limiting  drag  source.  This  is  an  important  issue  for  potential  methods, 
which  use  an  isentropic  flow  model.  It  has  been  shown  [76-78]  that  potential  methods  can 
be  simply  corrected  to  account  for  this  entropy  rise.  One  of  the  simplest  such  corrections 
has  been  applied  to  the  FPR  code  by  Bridgeman  et  al.  [37].  This  consists  of  applying  the 
nonisentropic  density  expression 


P  =  Pie 


-AS/R 


(56) 


where  pi  is  the  isentropic  density  given  by  Eq.  (35)  and  AS  is  the  shock-generated  entropy 
computed  by  the  Rankine-Hugoniot  relations.  It  can  be  shown  that  it  suffices  to  use  this 
expression  only  at  shock  points.  Equation  (56)  is  strictly  valid  only  for  steady  flows,  but  it 
has  produced  good  results  for  both  steady  and  unsteady  solutions.  In  general,  it  requires  an 
unusually  strong  shock  (for  rotors)  for  these  corrections  to  be  necessary.  Such  flows  often 
involve  shock-induced  separation  and  probably  require  at  least  a  TNS  flow  model.  Although 
the  use  of  a  strongly  interacting  boundary  layer  is  also  possible  these  methods  tend  not  to  be 
very  robust  for  strong  shocks.  Of  course,  all  these  analyses  encounter  the  turbulence  modeling 
problem.  Our  reliance  on  testing  will  not  soon  disappear  for  flows  with  strong  shocks,  but 
we  now  have  the  computational  tools  and  the  option  to  design  configurations  that  avoid  such 
situations. 

Figure  14  showed  a  comparison  of  the  computed  and  measured  torque  for  a  nonlifting 
rotor.  This  FPR  computation  employed  a  2-D,  noninteracting  integral  boundary-layer  model 
(Nash-MacDonald).  Subsequent  computations  [79,  80]  have  also  used  the  3-D  finite-difference 
boundary-layer  code  BL3D  [81].  For  this  simple  rectangular  rotor,  drag  results  are  very 
similar  for  the  two  boundary-layer  models,  the  integral  method  predicting  somewhat  lower 
values.  However,  beyond  drag  divergence  the  difference  is  inconsequential  because  most  of  the 
drag  is  due  to  inviscid  effects.  In  order  to  further  test  our  drag  prediction  ability,  nonlifting 
comparisons  have  also  been  performed  for  nonrectangular  rotors.  Figure  49  shows  a  sketch 
of  the  tested  rotor  and  two  pressure  distributions  plotted  as  functions  of  chord  and  profile 
height  (as  required  for  drag  integration).  It  can  be  seen  in  Fig.  50  that  the  computed  inviscid, 
spanwise  drag  distributions  for  rectangular  and  swept  planforms  are  very  different.  For  the 
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rectangular  case  the  drag  is  very  close  to  zero  except  near  the  tip.  Much  of  this  small  drag  is 
numerical.  However,  the  inviscid  drag  for  even  a  nonlifting  rotor  is  nonzero  due  to  acoustic 
radiation.  Figure  51  shows  a  comparison  of  computed  and  measured  torque  for  the  swept 
rotor,  using  the  integral  and  finite-difference  boundary-layer  models.  It  appears  that  fairly 
simple  boundary-layer  corrections  will  be  applicable  for  a  wide  range  of  conditions.  The  ability 
to  directly  integrate  surface  pressure  to  obtain  drag  indicates  that  we  should  also  be  able  to 
find  induced  drag  by  similar  means.  However,  this  depends  on  our  ability  to  accurately  predict 
the  wake. 


VI.  CONCLUDING  REMARKS 

Rotorcraft  aerodynamics  is  especially  rich  in  unsolved  problems,  and  for  this  reason  the 
need  for  interdependent  computational  and  experimental  studies  is  great.  Rotor  CFD  is  unique 
in  that  its  developers  have  frequently  been  experimentalists  as  well.  This  has  maintained  a 
needed  balance  between  our  ability  to  compute  and  our  ability  to  see  the  whole  problem. 
Considerable  progress  has  been  made  and  we  can  begin  to  think  of  rotor  CFD  as  a  viable 
engineering  tool  as  well  as  a  means  for  basic  studies. 

3-D  unsteady,  nonlinear  potential  methods  are  becoming  fast  enough  to  enable  their  use 
in  parametric  design  studies.  At  present,  combined  CAMRAD/FPR  analyses  for  a  complete 
trimmed  rotor  solution  can  be  performed  in  about  an  hour  on  a  CRAY  Y-MP  (or  ten  minutes, 
with  multiple  processors).  These  computational  speeds  indicate  that  in  the  near  future  many 
of  our  large  CFD  problems  will  no  longer  require  a  supercomputer.  It  is  also  becoming  clear 
that  potential-based  methods  are  more  capable  than  we  had  previously  supposed.  The  ability 
to  convect  circulation  is  routine  for  integral  methods,  but  only  recently  have  we  discovered  how 
to  do  the  same  with  differential  methods.  With  the  HELIX-I  code  it  is  possible  to  compute  an 
entire  hover  performance  polar  (about  five  full  3-D,  supercritical  flow  computations,  including 
boundary-layer  and  free-wake)  in  about  six  hours.  Steady,  viscous  airfoil  computations  (for 
example,  with  ARC2D)  for  flows  with  no  major  separation  can  be  done  in  minutes.  These 
could  greatly  augment  our  still  much-used  airfoil  data  base. 

It  is  clear,  then,  that  the  differential  CFD  rotor  analyses  are  poised  to  enter  the  engineer¬ 
ing  workplace.  Integral  methods  already  constitute  a  mainstay.  Although  much  development 
is  still  required,  the  major  need  now  is  for  skillful  users  who  can  apply  these  tools  to  their 
own  individual  problems.  Ultimately,  it  is  these  users  who  will  integrate  CFD  into  the  entire 
engineering  process  and  provide  a  new  measure  of  confidence  in  design  and  analysis. 

It  should  be  recognized  that  the  above  classes  of  analyses  do  not  include  several  major 
limiting  phenomena  (especially  dynamic  stall),  which  will  continue  to  require  empirical  treat¬ 
ment  because  of  computational  time  constraints  and  our  limited  physical  understanding.  Such 
empirical  treatment  should  be  included,  however,  into  our  developing  CFD,  engineering-level 
analyses.  Moreover,  it  is  probable  that  in  the  near  future  CFD  will  be  reliable  enough  to  pro¬ 
vide  a  new  source  of  empirical  information  with  which  to  supplement  physical  measurements. 
We  can  expect  to  be  able  to  visualize  effects  and  test  ideas  in  ways  that  are  not  possible  with 
physical  testing.  It  is  likely  that  properly  constructed  flow  models  containing  corrections  from 
physical  testing  will  be  able  to  fill  in  unavoidable  gaps  in  our  experimental  data  base,  both  for 
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basic  studies  and  for  specific  configuration  testing.  For  these  kinds  of  applications,  computa¬ 
tional  cost  is  not  an  issue.  For  rotorcraft,  testing  and  computation  will  become  increasingly 
and  truly  interdependent,  and  the  extent  of  this  integration  should  be  an  important  measure 
of  their  effectiveness. 

Finally,  we  should  recognize  that  although  rotorcraft  are  probably  the  most  complex 
of  aircraft,  the  rotorcraft  engineering  community  is  very  small  compared  to  the  fixed-wing 
community.  Likewise,  rotorcraft  CFD  resources  can  never  achieve  fixed-wing  proportions  and 
must  be  used  wisely.  Therefore  we  must  glean  the  fixed-wing  work  for  many  of  our  basic 
methods.  This  approach  has  its  limits,  though,  because  rotor  needs  are  unique  and  cannot  be 
met  without  much  original  thinking.  This  is  a  fertile  field  with  much  yet  to  accomplish. 
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Fig.  2  Panel  computation  of  a  complex  rotor/wake/fuselage  flow  using  VSAERO.  Details  of 
the  overall  wake  structure  of  a  V-22  in  a  typical  transition  flight  condition.  Airspeed  45  knots, 
a  =  -1.5°,  Ct  =  0.015  [7]. 
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Fig.  4  A  local  model  for  self-induction  of  a  vortex  element. 
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Fig.  5  The  free-wake  performance  prediction  of  an  S-76  rotor  [17]. 


Fig.  8  Grid  and  boundary  conditions  for  a  local  rotor-blade  computation. 
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Fig.  9  A  model  problem,  (a)  Solving  for  the  plunging-pitching  motion  of  an  airfoil  using  the 
low-frequency,  transonic,  small-perturbation  equation,  (b)  Difference  approximations  required 
for  the  model  problem. 
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Fig.  10  The  transonic  flow  on  a  rotor  blade,  measured  and  computed  using  the  low-frequency 
transonic  small-perturbation  equation,  showing  the  effect  of  transonic  unsteadiness  due  to  the 
azimuthally  varying  Mach  number  [33] . 


Fig.  11  A  comparison  of  unsteady  and  quasi-steady,  3-D  transonic  rotor  computations  using 
the  TFAR2  nonconservative,  full-potential  code  [27]. 
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Fig.  12  A  novel  method  for  obtaining  high-speed  rotor  flow-field  data  for  code  validation: 
interferograms  of  a  transonic  rotor  at  various  azimuth  angles  [36] . 
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Fig.  13  A  comparison  of  measured  and  computed  near-surface  pressure  coefficients.  Data 
obtained  from  tomographic  reconstruction  of  interferograms  [36]. 


Fig.  14  Measured  and  computed  torque  for  a  nonlifting  rotor. 
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Fig.  15  A  finite-difference  model  embedded  in  a  global  hover  fiow,  showing  the  disparity 
between  a  typical  blade-oriented  grid  and  the  total  hovering  flow  field. 
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Fig.  16  Near- field  blade/vortex  modeling  methods  for  use  with  hybrid  methods. 


Mj  =  0.877, 0- =  8”,  AR  =  6.0  -  Conservative  full  potential  [30] 

^  -  Non-conservallve  full  potential  [41] 

•  Experiment  [42]  .  Small  disturbance  [38] 


Fig.  17  Measured  and  computed  flow  on  a  hovering  rotor — a  comparison  of  various  hybrid 
potential  methods  for  predicting  high-speed,  hovering  flows. 
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Fig.  18 


A  hybrid  rotor  flow  prediction  scheme  for  forward  flight,  using  alpha  coupling  [43]. 


Fig.  19  A  flow  diagram  for  the  coupled  CAMRAD  and  FPR  codes  [43]. 
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.  20  A  comparison  of  surface  pressure  data  from  the  ONERA  3-blade  model  rotor  with 
MRAD/FPR  results  [44]. 


52 


Fig.  22  A  comparison  of  hover  surface  pressure  data  with  results  from  a  unified  rotor/wake 
computation  using  a  centered  Euler  scheme  [51]. 


Fig.  23  The  dissipation  of  vorticity  in  a  unified  rotor/wake  computation  using  an  upwind 
Euler  scheme  [53]. 
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Fig.  24  Alternate  models  of  the  shed  wake  showing  a  wake  reformulation  that  permits  free 
circulation  convection  in  a  potential  method  [55]. 


Fig.  25  A  computed  wake  structure  for  a  hovering  rotor  using  the  HELIX-I  vorticity- 
convecting,  full-potential  code  [56]. 
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Fig.  28  A  comparison  of  computed  and  measured  rotor  hover  performance  for  a  4-blade  rotor 
[56]. 
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Fig.  29  A  comparison  of  computed  and  measured  rotor  hover  performance  for  the  Boeing 
Model  360  rotor  [58]. 


Fig.  30  The  computed  free  convecting  wake  for  the  Boeing  model  360  rotor  [58]. 
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Fig.  32  Lift  and  moment  variations  for  a  dynamic  stall  using  an  integral-diflferential  scheme 
[62]. 
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Fig.  35  The  lift  variation  for  a  dynamically  stalling  airfoil  with  constant  pitch  rate:  the  effect 
of  pitch  rate  [66]. 
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Fig.  36  The  lift  variation  for  a  dynamically  stalling  airfoil  with  constant  pitch  rate:  the  effect 
of  Mach  number  [66]. 
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Fig.  37  A  particle  trajectory  visualization  of  the  steady  stall  of  a  nonrectangular  rotor  blade, 
computed  using  TNS  [67]. 


(a)  Round  Up 


(b)  Square  tip 


Fig.  38  The  formation  and  liftoff  of  the  tip  vortex  for  a  rectangular  wing.  Mo©  =  0.17, 
a  =  11.8°,  and  Re  =  2  x  10®  [70]. 


(b) 


Fig.  39  Tip>-vortex  studies:  a  comparison  of  computed  and  experimental  surface  oil  flow  pat¬ 
terns.  Moo  =  0.17,  Re  =  2  X  10®  [70].  (a)  Round  tip-experiment,  (b)  Round  tip-calculations. 
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Fig.  40  Surface  pressure  distributions  in  the  tip  region  of  a  rectangular  wing  with  different 
tip  caps.  Moo  =  0-17,  q  =  11.8°,  and  Re  =  2  x  10®  [70]. 
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Fig.  41  Instantaneous  surface-pressure  distributions,  adaptive  grid,  and  pressure  contours 
during  an  airfoil  interaction  with  a  convecting  vortex.  NACA0012,  Moo  =  0.8,  a  =  0°, 
Re  =  2  X  10®,  and  closest  approach  distance  =  -0.26  chords  [71]. 
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Fig.  42  An  Euler/TNS  grid  for  use  in  airfoil/vortex  interaction  studies  [73]. 


Fig.  43  Vorticity  contours  for  a  head-on  blade/vortex  interaction  (BVI)  computed  using  a 
high-order  upwind  TNS  scheme.  Moo  =  0.536  [73]. 
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Fig.  45  Pressure  contours  for  a  supercritical  BVI  [74]. 


Fig.  46  A  rotor  BVI  test  setup  for  the  acquisition  of  surface  pressure  data 


Fig.  47  A  comparison  of  measured  BVI  surface  pressures  with  computed  results  from  a 
high-order  upwind  TNS  scheme  [73]. 
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Fig.  48  A  comparison  of  measured  BVI  surface  pressures  with  3-D,  full-potential  results  for 
a  supercritical  oblique  interaction  [47] . 


Fig.  49  Computed  surface  pressures  for  a  swept  rotor  at  two  radii,  showing  the  pressure- 
thickness  plot  for  drag  integration  [81]. 


Fig.  50  The  computed  spanwise  drag  distribution  for  swept  and  rectangular  blades  [81]. 


Tip  mach  number 


Fig.  51  Computed  and  measmred  torque  for  a  swept-tip  rotor  in  hover  [81]. 


72 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Pubiic  reporting  burdbn  for  this  cofloction  of  information  is  astimatsd  to  avarags  t  hour  par  rasponsa.  induding  tha  tima  for  rayiawing  instructions,  saarching  axisting  data  sourcas, 
gathering  and  maintaining  tha  data  needed,  and  completing  and  reviewing  tha  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information.  Including  suggestions  tor  reducing  this  burden,  to  Washington  Headqu^ers  Services.  Directorate  for  information  Operations  and  Reports.  12tS  Jefferson 
Davis  Highway.  Suite  1204.  Arlington.  VA  22202-4302.  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduction  Project  (0704-0188).  Washington.  DC  20503. 

1.  AGENCY  USE  ONLY  fLajv*  b/an*;  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

March  1992  Technical  Memorandum 

4.  TITLE  AND  SUBTITLE 

The  Application  of  CFD  to  Rotary  Wing  Flow  Problems 

5.  FUNDING  NUMBERS 

505-61-51 

6.  AUTHOR(S) 

F.  X.  Caradonna 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  AOORESS(ES) 

Ames  Research  Center,  Moffett  Field,  CA  94035-1000  and 
Aeroflightdynamics  Directorate,  U.S.  Army  Aviation  Systems 
Command,  Ames  Research  Center,  Moffett  Field,  CA  94035-1099 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

A-90110 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

National  Aeronautics  and  Space  Administration 

Washington,  D.  C.  20546-0001  and  U.S.  Army  Aviation  Systems 
Command,  St.  Louis,  MO  63120-1798 

10.  SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

NASA  TM- 102803 
USAAVSCOM  TR-92-A-004 

11.  SUPPLEMENTARY  NOTES 

Point  of  Contact:  F.  X.  Caradonna,  Ames  Research  Center,  MS  215-1,  Moffett  Field,  CA  94035-1000 
(415)  604-5902  or  FTS  464-5902 

This  paper  was  originally  prepared  as  lecture  notes  for  an  Agard  course,  “Aerodynamics  of  Rotoicraft.” 

12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Unclassified-Unlimited 

Subject  Category  -  02 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  rMav/mum  200  worcfs; 

This  report  is  a  review  intended  to  serve  as  source  and  reference  material  for  the  CFD  phase  of  a  short  course  on 
rotorcraft  aerodynamics.  The  review  is  mainly  oriented  toward  engineering  application  methods.  Therefore  emphasis 
is  placed  on  potential  metliods  and  their  various  types  and  uses.  However,  the  application  of  viscous  codes  is  also 
discussed.  These  methods  are  discussed  in  the  context  of  the  various  rotorcraft’-specific  problems — including  wake 
prediction,  bladeAortex  interactions,  transonic  problems,  dynamic  stall  and  other  separation  problems. 

The  exposition  includes  a  discussion  of  the  various  flow  equations  and  the  physical  approximations,  which  they 
embody.  Basic  methods  of  solving  these  are  presented — especially  integral  and  difference  approaches  to  the  various 
potential  equations.  The  treatment  of  total  rotor/wake  problems  by  hybrid  difference/integral  schemes  are  discussed. 
These  include  both  hover  flow  and  forward  flight  methods  requiring  a  consideration  of  trim.  Unified  schemes,  which 
solve  an  entire  rotor/wake  flow  on  a  single  grid,  are  also  treated.  One  of  these,  a  vorticity  convecting  full-potential 
approach,  is  capable  of  predicting  complete  hover  performance  polars.  Many  comparisons  are  made  with  data  from 
flight  and  model  tests,  to  demonstrate  the  efficacy  of  the  different  flow  treatments.  Finally,  advanced  flow  topics, 
requiring  full  viscous  solutions,  are  shown  in  order  to  demonstrate  the  future  possibilities  and  opportunities  for  CFD 
as  newer  methods  mature  and  become  practical. 

14.  SUBJECT  TERMS 

Rotorcraft  aerodynamics.  Wakes,  CFD,  Transonic  flows.  Viscous  flows 

IS.  NUMBER  OF  PAGES 

74 

16.  PRICE  CODE 

A04 

17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION 
OF  REPORT  OF  THIS  PAGE 

Unclassified  Unclassified 

19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

20.  LIMITATION  OF  ABSTRACT 

NSN  7540-01 -280'5500  Standard  Form  298  (Rav.  2'89) 


Pr«scrlb#tf  by  ANSI  Sttf.  Z9S-10 


2M-102 


